Kf1 - Tillampad Matematik, TMA 225 - 2007

Studio 3: FEM i 1D.

Börja med att läsa igenom det utdelade materialet om Robinrandvillkor.

Dagens uppgift är att arbeta med och utveckla programmet PoissonSolver1D till en komplett FEM-lösare för stationära (d.v.s. icke-tidsberoende) reaktions-diffusionsproblem med Robinrandvillkor i en rumsdimension. Eftersom kunskaperna från den förra studioövningen om beräkning av L2-projektion är till stor nytta för att förstå hur FEM-lösaren fungerar, är det viktigt att du tagit till dig detta material innan du börjar med dagens övning.

Huvudmomenten i PoissonSolver1D är: 

  - Skapa en partition av [x_min, x_max] i N-1 stycken delintervall.

  - Assemblera matriser (S, M, R) och vektorer (v, r).

  - Lös det linjära ekvationssystemet.

  - Visualisering, dvs plotta lösningen.
Dessa steg är implementerade på följande sätt:

Först skapas en partition av [x_min, x_max] i huvudprogrammet PoissonSolver1D, vilket därefter anropar funktionen PoissonAssembler1D, i vilken assembleringen av matriser S, M och R och vektorer v och r görs.

Genom att låta
>> A = S + M + R;
och
>> l = v + r;
fås lösningsvektorn U som
>> U = A\l;
Lösning av det linjära ekvationssystemet och plottning av lösningen görs i PoissonSolver1D.

Gå igenom koden i PoissonSolver1D.m noggrant, så du förstår programmets struktur. Denna kod behöver du inte göra något med själv.

Fortsätt med att gå igenom koden i PoissonAssembler1D.m. Assembleringen av lastvektorn, v, är identisk med hur vi gjorde när vi beräknade högerledsvektorn i fallet med L2-projektion, och finns därför redan inlagd. Också assembleringen av massmatrisen, M, är (så när som på koefficienten c) densamma som för L2-projektion, och är även den inlagd.

Det du alltså själv måste göra är att, på de platser i koden där det står ???, komplettera PoissonAssembler1D.m så att styvhetsmatrisen, S, assembleras, och randbidragen till ekvationssystemet, matrisen R och vektorn r, beräknas. (Jämför med hur assembleringen av M och b görs: För beräkning av M används Simpsons formel och för beräkning av b används Trapetsregeln. Du kan förslagsvis välja Trapetsregeln för beräkning av styvhetsmatrisen. Jämför också med det utdelade materialet om Robinrandvillkor.)

Prova nu att lösa exemplet som finns beskrivet i PoissonSolver1D.m Kör programmet med kommandot >>PoissonSolver1D. Det är mycket viktigt att testa ditt program genom att jämföra din beräknade lösning med den exakta lösningen (vilken är känd i detta fall), på det sätt som finns beskrivet i slutet av PoissonSolver1D.m.

Lycka till!!!

Tillbaka till kurshemsidan.


Editor: Nils Svanstedt
Last modified: 2007-03-30