K2/Bt2 - TMK: Tillämpad matematik K (TMA 682) - 2003/2004

Datorstudioövning 2: L2-projektion och kvadratur.

På dagens studioövning börjar vi med att studera några olika kvadraturregler, som vi kommer att använda för numerisk beräkning av integraler. Dessa regler känner du säkert till sedan tidigare, t.ex. från ALA-B. Vi skall tillämpa dessa regler för numerisk beräkning av L2-projektion. Detta är en metod för att approximera funktioner med styckvis kontinuerliga polynom.

Starta med att hämta de två filerna LoadVector.m och MassMatrix.m från kursens hemsida. Lägg dessa filer i din hemkatalog och starta Matlab.
För att beräkna L2-projektionen, &Pi f, till en given funktion f(x) på intervallet [a,b]. Skriv därför en funktionsfil, f.m, som definierar f(x) t.ex: 


  function y = f(x)
  y = x.^2;
Prova att plotta f(x) på intervallet [0,1]: 

  >> fplot('f', [0 1])
För att beräkna &Pi f måste vi först beräkna massmatrisen, M, och lastvektorn, v. Här använder vi oss av en styckvis linjär approximation med bas bestående av hattfunktioner.
Innan du anropar LoadVector som beräknar högerledet, måste du skapa en vektor, som innehåller nodpunkterna i partitionen av [a,b]. T.ex. för att beräkna högerledet på en likformig indelning av [0,1] i fyra delintervall skriver du: 

  >> p = 0:1/4:1;
  >> v = LoadVector(p)
Beräkna därefter massmatrisen genom att anropa: 

  >> M = MassMatrix(p)
Hur ser M ut? Jämför med de analytiska uttrycken som du kan beräkna för hand.

Som du kanske märker överensstämmer inte M med den analytiskt beräknade massmatrisen. Det beror på att kvadraturregeln, trapetsregeln, som används i koden bara är exakt för polynom av grad mindre än eller lika med 1. Integranderna i massmatrisen är dock polynom av grad 2, dvs. med trapetsregeln får vi bara en approximativ massmatris.

Nu kan du beräkna &Pi f:s komponentvektor, c, genom att lösa det linjära ekvationssystemet M c = v, enligt: 


  >> c = M\v
Plotta f(x) och &Pi f i samma figur: 

  >> figure                   
  >> a = p(1);                % a = x_0 = första elementet i p
  >> b = p(end);              % b = x_N =  sista elementet i p
  >> fplot('f', [a b])   
  >> hold on                  % lås figuren
  >> plot(p, c, 'rd-')        % plotta &Pi f 'v'.
Hmm... ser inte detta ut som den vanliga nodinterpolanten till f(x)? Prova att plotta också denna i samma figur: 

  >> plot(p, f(p), 'go-')
Hur kommer detta sig? Jo, det beror på att vi med trapetsregeln erhåller en approximativ massmatris.

Extra uppgift:

Förbättra approximationen av &Pi f genom att använda en noggrannare kvadratur, t.ex. Simpsons formel. Det viktigaste är att byta kvadraturregel i funktionen MassMatrix, eftersom integranden där har en stor andraderivata.

Jämför ånyo med de analytiska uttrycken. Stämmer det nu?
Upprepa slutligen jämförelsen ovan mellan L2-projektionen och nodinterpolanten genom att rita dem i samma figur och studera skillnaden. Det är lättast att se om du väljer ett litet antal delintervall. Tänk efter hur de är definierade och hur detta påverkar deras utseende.

Lycka till!

Tillbaka till kurshemsidan.


Editor: Fredrik Bengzon
Last modified: 2003-09-07