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

Datorstudioövning 3: Styvhetsmatris. Robinrandvillkor. Börja med att läsa igenom det utdelade materialet om Robinrandvillkor. Dagens uppgift är att bygga ut programmet MyFirstPoissonSolver 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 MyFirstPoissonSolver är: - Skapa en partition av [x_min, x_max] i N-1 stycken delintervall. - Assemblera matriser (A, Mc, R) och vektorer (b, rv). - Lös det linjära ekvationssystemet. - Plotta lösningen. Dessa steg är implementerade på följande sätt: Först skapas en partition av [x_min, x_max] i huvudprogrammet MyFirstPoissonSolver, vilket därefter anropar funktionen MyFirstPoissonAssembler, i vilken assembleringen av matriser och vektorer görs. MyFirstPoissonAssembler anropar i sin tur funktionerna a, c, f, vilka specificerar koefficienterna i differentialekvationen, samt gamma, g_D, g_N, vilka specificerar Robinrandvillkoren. Lösning av det linjära ekvationssystemet och plottning av lösningen görs slutligen i MyFirstPoissonSolver. Du öppnar huvudprogramsfilen MyFirstPoissonSolver.m i Matlabs editor genom att skriva >> edit K:\Matkfkb\0102\detb-a-k-kf-1\MyFirstPoissonSolver Spar programmet i ett av dina egna bibliotek. Fortsätt nu med att öppna, samt spara hos dig själv, funktionsfilerna MyFirstPoissonAssembler.m, a.m, c.m, f.m, gamma.m, g_D.m, och g_N.m, vilka samtliga ligger i kursbiblioteket K:\Matkfkb\0102\detb-a-k-kf-1\. Gå igenom koden i MyFirstPoissonSolver.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 MyFirstPoissonAssembler.m. Assembleringen av lastvektorn, b, är identisk med hur vi gjorde när vi beräknade högerledsvektorn, b_vek, i fallet med L2-projektion, och finns därför redan inlagd. Också assembleringen av massmatrisen, Mc, ä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 MyFirstPoissonAssembler.m så att styvhetsmatrisen, A, assembleras, och randbidragen till ekvationssystemet, matrisen R och vektorn rv, beräknas. (Jämför med hur assembleringen av Mc och b görs: För beräkning av Mc 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 MyFirstPoissonSolver.m: Kontrollera först att filerna a.m, c.m, f.m, gamma.m, g_D.m, och g_N.m är korrekt definierade, och kör slutligen programmet med kommandot >> MyFirstPoissonSolver. 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 MyFirstPoissonSolver.m. Lycka till!!! Tillbaka till kurshemsidan.

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