Kf1 - Differentialekvationer och tekniska beräkningar del A, TMA 225 - 2003

Datorstudioövning 2: Kvadratur. L2-projektion.

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. En del av dessa regler känner du säkert till sedan tidigare, t.ex. från ALA-B. Vi fortsätter sedan med att tillämpa dessa kunskaper på numerisk beräkning av L2-projektion.

Börja med att läsa igenom det utdelade materialet om kvadratur i en rumsdimension, och försök med dess hjälp lösa de första två problemen i veckans Räkneuppgifter.

Vi skall nu beräkna L2-projektionen, Phf, till en given funktion, f, i L2(a,b). Skapa och spara först en funktionsfil, f.m, som definierar f. Den kan t.ex. se ut så här:


  function y = f(x)
  y = x.^2;   % Jämför gärna med Problem 3. på veckans Räkneuppgifter.

Det är bäst att använda point-wise operationer som ovan, för att t.ex. kunna plotta f med kommandot fplot. Prova att plotta f på intervallet [0,1]:

  >> fplot('f', [0 1])

För att beräkna Phf måste vi först beräkna massmatrisen, M, och högerledsvektorn, b_vek. (Vi kallar högerledsvektorn för b_vek för att undvika sammanblandning med intervallets högra ändpunkt, som vi ju ofta kallar b!)

Vi börjar med b_vek. Skriv


  >> edit K:\Matkfkb\0102\detb-a-k-kf-1\MinHLAssemblerare

så öppnas en funktion som beräknar b_vek genom att approximera integralerna med trapetsregeln. Spara MinHLAssemblerare.m i ett av dina egna bibliotek (du kommer att få modifiera den lite senare).

Gå mycket noga igenom koden och försök förstå hur den fungerar: Tänk på att syftet är att beräkna vektorn b_vek, vilken har N+1 komponenter, och där varje komponent fås som en integral av f gånger motsvarande hatt-funktion. Funktionen skall alltså beräkna N+1 stycken integraler! Eftersom hatt-funktionerna endast är skilda från noll på två delintervall ("hatten" längst till vänster och "hatten" längst till höger endast på ett delintervall!) får dock varje integral endast bidrag från (högst) två delintervall. Ett sätt att beräkna b_vek vore att räkna ut ett element (d.v.s. en integral) i taget genom att addera bidragen från de intervall där integranden är skild från noll. I MinHLAssemblerare.m är dock beräkningen organiserad på ett annat sätt: Funktionen loopar över de N delintervallen i tur och ordning, och räknar ut samtliga integralbidrag från ett delintervall i taget. Ett fixt delintervall bidrar till exakt två av integralerna, eftersom alla hatt-funktioner utom två är identiskt lika med noll på intervallet.

Denna teknik, att addera bidragen till vektorns element från ett delintervall i taget, kallas för assemblering.

Innan du anropar MinHLAssemblerare måste du skapa en radvektor som innehåller nodpunkterna i partitionen av [a,b]. För att skapa en likformig indelning av [0,1] i fyra stycken delintervall skriver du:


  >> p = 0:1/4:1;

Nu kan du beräkna högerledsvektorn genom att anropa MinHLAssemblerare:

  >> b_vek = MinHLAssemblerare(p)

Nästa steg är att beräkna massmatrisen: På samma sätt som ovan öppnar du och sparar i ett eget bibliotek funktionsfilen MinMassmatrisAssemblerare.m (som också den ligger i biblioteket K:\Matkfkb\0102\detb-a-k-kf-1\). Gå noga igenom också denna kod. Även MinMassmatrisAssemblerare assemblerar matriselementen genom att loopa över delintervallen i tur och ordning, och använder trapetsregeln för att approximera integralerna på varje delintervall. Notera att det nu är (N+1)^2 integraler som skall beräknas! Ett fixt delintervall bidrar till exakt fyra av matriselementen (integralerna), eftersom alla hatt-funktioner utom två är identiskt lika med noll på intervallet.

Beräkna nu massmatrisen genom att anropa MinMassmatrisAssemblerare:


  >> M = MinMassmatrisAssemblerare(p)

Hur ser M ut? Jämför med dina analytiska uttryck i Problem 4. i denna veckas Räkneuppgifter!

Som du kanske märker överensstämmer inte M med den analytiskt beräknade massmatrisen. Det beror på att trapetsregeln bara är exakt för polynom av grad mindre än eller lika med 1, men integranderna i massmatrisen är polynom av grad 2. Så med trapetsregeln får vi bara en approximativ massmatris. (Det vi får är faktiskt en diagonal matris som kallas för den "lumpade" massmatrisen!) Lite senare skall vi använda en noggrannare kvadraturregel (t.ex. Simpsons formel) som ger exakt resultat i detta fall, men tills vidare nöjer vi oss med den approximativa massmatrisen.

Nu kan du beräkna Phf:s komponentvektor, c, genom att lösa det kvadratiska, linjära ekvationssystemet M c = b_vek:


  >> c = M\b_vek

Plotta f och Phf i samma figur:

  >> figure                   % öppna ett nytt figurfönster
  >> a = p(1);                % a = x_0 = första elementet i p
  >> b = p(length(p));        % b = x_N =  sista elementet i p
  >> fplot('f', [a b], 'k')   % plotta f över [a,b] ('k' => svart färg)
  >> hold on                  % "lås" figuren
  >> plot(p, c, 'bd-')        % plotta Phf ('b' => blå färg,
                              %    'd' => rita "diamanter" i datapunkterna,
                              %    '-' => sammanbind datapunkterna med
                              %           heldragna linjer)

Hmm... Ser inte detta ut som den vanliga nod-interpolanten till f?! Prova att plotta också denna i samma figur:

  >> plot(p, f(p), 'ro-')     % 'r' => röd färg

Hur kommer detta sig? Jo, vi observerade tidigare att vi med trapetsregeln erhåller en approximativ massmatris, nämligen den "lumpade" massmatrisen, som är en diagonal matris med följande element på huvuddiagonalen: [h/2 h h h ... h h h h/2]. (Om vi antar att indelningen är likformig med delintervall av längd h.) Kommentar: Vi kan i detta sammanhang nämna att om vi istället använt mittpunktsformeln som kvadraturregel, och haft en likformig indelning, vi erhållit en approximativ massmatris som ej är inverterbar. Mittpunktsformeln är alltså ej en lämplig kvadraturregel i detta fall.

När vi vidare använder trapetsregeln för att approximera b_vek blir resultatet b_vek = [h/2*f(x_0); h*f(x_1); h*f(x_2); ...; h*f(x_N-2); h*f(x_N-1); h/2*f(x_N);]. Jämför med Problem 1. (e) i veckans Räkneuppgifter.

Lösningen till ekvationssystemet M c = b_vek blir därför i detta fall inget annat än c = [f(x_0); f(x_1); f(x_2); ...; f(x_N-2); f(x_N-1); f(x_N);], d.v.s., den approximation av L2-projektionen, Phf, som man får då man approximerar integralerna med trapetsregeln är nodinterpolanten! Kommentar: Notera hur enkelt det linjära ekvationssystemet är att lösa eftersom den "lumpade" massmatrisen är diagonal: Detta är en fördel med denna approximation.

För att förbättra approximationen av Phf får du använda en noggrannare kvadraturregel, antingen Simpsons formel eller Gauss-Legendres 2-punktsformel. Det viktigaste är att byta kvadraturregel i funktionen MinMassmatrisAssemblerare, eftersom integranden där har en stor andraderivata, proportionell mot 1/h^2, medan dess tredjederivata är noll. Vi vinner därmed mycket i noggrannhet genom att byta till en högre ordningens kvadraturmetod. Vi erhåller faktiskt exakta värden på massmatriselementen både med Simpsons formel och Gauss-Legendres 2-punktsformel, eftersom integranden är ett andragradspolynom.

Försök nu alltså att byta till en noggrannare kvadraturregel genom att modifiera MinMassmatrisAssemblerare (och om du vill förstås också MinHLAssemblerare). Jämför ånyo med de analytiska uttrycken i Problem 4. i denna veckas Räkneuppgifter! 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.

Tips: Om du även modifierat MinHLAssemblerare kan du använda denna till att kontrollera högerledsvektorerna du skall beräkna för hand i Problem 3. (a) - (b) i veckans Räkneuppgifter. Du använder då förstås >> p = 0:1:1; respektive >> p = 0:1/2:1;. (Eftersom f(x) är ett andragradspolynom blir integranden ett tredjegradspolynom, vilket både Simpsons formel och Gauss-Legendres 2-punktsformel integrerar exakt!)

Lycka till!

Tillbaka till kurshemsidan.


Editor: Georgios Foufas
Last modified: 2003-03-24