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

Datorstudioövning 5: Tidsberoende problem. På dagens studioövning skall vi studera tidsberoende reaktions-diffusionsproblem, d.v.s., vi tänker oss nu att lösningen u = u(x, t) kan bero på både rum x och tid t. Framställningen nedan är i första hand avsedd som ett stöd för Matlab-implementering, och är inte fullständig. Se dina föreläsningsanteckningar samt det utdelade materialet om Tidsberoende problem, för ytterligare detaljer rörande diskretisering i rum och tid. Modellering. För att få en fysikalisk känsla för ekvationen vi skall lösa kan vi tänka oss att u är koncentrationen hos ett reagerande ämne, och börjar med att studera det stationära fallet som vi tidigare betraktat. I detta fall måste i varje ``punkt'' en lika stor mängd av ämnet tillföras som bortföras. Annars skulle ju koncentrationen inte vara konstant i tiden. Betrakta ett litet intervall [x, x + dx] med längd dx. Nettotillförseln per tidsenhet av ämnet är a*u'(x + dx) - a*u'(x) - c*u(x)*dx + f(x)*dx. Den första termen beskriver den mängd som diffunderar genom intervallets högra ändpunkt och den andra termen den mängd som diffunderar genom intervallets vänstra ändpunkt. (Tänk på att -a*u' är positivt för ett flöde riktat åt höger.) Den tredje termen beskriver den mängd som reagerar och den fjärde termen beskriver produktion av ämnet. I ett stationärt tillstånd måste nettotillförseln vara lika med noll, vilket efter omflyttning av termer och division med dx ger: -(a*u'(x + dx) - a*u'(x))/dx + c*u(x) = f(x). Låter vi slutligen dx gå mot noll erhåller vi det stationära reaktions-diffusionsproblemet -(a*u')' + c*u = f, som vi ju känner sedan tidigare. Antag nu att det ej råder balans mellan bortförd och tillförd mängd av ämnet i en punkt x vid tiden t. Efter en (kort) tid dt har det då skett en nettoökning u(x, t + dt) - u(x, t) ~ ((a*u')' - c*u + f)*dt av koncentrationen. Dividerar vi med dt och låter dt gå mot noll erhåller vi det tidsberoende reaktions-diffusionsproblemet du/dt - (a*u')' + c*u = f, som vi nu skall studera. Detta kallas för en partiell differentialekvation (PDE) eftersom det förekommer derivator m.a.p. flera oberoende variabler: I detta fall en rumsvariabel och en tidsvariabel. Liksom tidigare betraktar vi Robin-randvillkor, vilka modellerar att det diffusiva flödet, -a*u', genom randpunkterna är proportionellt mot skillnaden, u - g_D, mellan koncentrationen precis vid randen och omgivningens koncentration (g_D). Proportionalitetskonstanten gamma är ett mått på randens permeabilitet. (Vi kan också tänka oss att vi på något sätt ``tvingar på'' ett flöde genom randen, vilket vi liksom tidigare modellerar med funktionen g_N.) Differentialekvation. Ekvationen vi skall lösa är alltså: du/dt - (a*u')' + c*u = f, x_min < x < x_max, 0 < t < T, a*u' = gamma*(u - g_D) + g_N, x = x_min, 0 < t < T, -a*u' = gamma*(u - g_D) + g_N, x = x_max, 0 < t < T, u(x, 0) = u0(x), x_min < x < x_max, t = 0. Kommentar: Att det diffusiva flödet i randvillkoret i punkten x = x_min har ombytt tecken beror på att vi har konventionen att ett utåtriktat flöde skall ha ett positivt tecken. Eftersom, som nämnts ovan, -a*u' är positivt för ett flöde riktat åt höger byter vi därför tecken i x = x_min. Vi önskar alltså bestämma lösningen u(x, t) i alla punkter x mellan x_min och x_max, för alla tider t fram till sluttiden T. Observera att eftersom detta är ett tidsberoende problem vi måste ange ett begynnelsevillkor, d.v.s. en funktion u0(x) som talar om vad lösningen är vid tiden t=0. I princip kan samtliga koefficienter i differentialekvationen, a, c, och f, och randvillkorsfunktioner, gamma, g_D och g_N, bero både på rum x och tid t, men vi nöjer oss tillsvidare med att betrakta fallet f = f(x, t), och antar att övriga data endast beror av x. Semi-diskretisering: Diskretisering i rummet. Vi börjar med att diskretisera i rummet. Vi söker alltså en approximation, U(x,t), till u(x, t) sådan att U(x,t) för varje fix tidpunkt t är en styckvis linjär, kontinuerlig funktion av x. För att åstadkomma detta gör vi följande ansats: U(x,t) = xi_1(t)*phi_1(x) + ... + xi_N(t)*phi_N(x) Skillnaden mot ansatsen i det stationära fallet är att vikterna, xi_j(t), nu är tidsberoende. Eftersom vi har Robin-randvillkor summerar vi över samtliga noder, också randnoderna. (Om man däremot har homogena Dirichlet-randvillkor, u = 0, tar man inte med basfunktionerna som hör till randnoderna.) Genom att variationsformulera differentialekvationen, och därefter ``byta'' u mot U samt välja hatt-funktionerna phi_i som testfunktioner, erhålls följande system av N stycken ordinära differentialekvationer (ODE): (Jfr. Problem 5 i veckans räkneuppgifter.) M*d(xi)/dt + (A + Mc + R)*xi(t) = b(t) + rv, 0 < t < T, xi_j(0) = u0(x_j), j = 1, ..., N. Ur detta system av ODE kan vikt-funktionerna xi(t) = [xi_1(t); xi_2(t); ...; xi_N(t)], vilka fullständigt bestämmer U(x,t), beräknas. Notera att begynnelsevillkoret xi_j(0) = u0(x_j) innebär att U(x,0) = u0(x_1)*phi_1(x) + ... + u0(x_N)*phi_N(x), d.v.s., vi väljer U(x,0) som nodinterpolanten till u(x, 0). Notera också att eftersom vi antagit att f = f(x, t) så kommer lastvektorn b = b(t) att bero på tiden. Vi har ju bara integrerat i rumsled ännu så länge! Kommentar: Det som skiljer fallet med Robin-randvillkor från fallet med homogena Dirichletrandvillkor, är förekomsten av termerna R*xi(t) samt rv. Matrisen R och vektorn rv är precis desamma som i det stationära fallet, och uppträder på ett identiskt vis. Diskretisering i tiden. Bakåt Euler. För att förenkla notationen inför vi nu beteckningarna A_tot = A + Mc + R samt b_tot(t) = b(t) + rv. Systemet av ODE kan då skrivas: M*d(xi)/dt + A_tot*xi(t) = b_tot(t), 0 < t < T, xi_j(0) = u0(x_j), j = 1, ..., N. Kommentar: Notera att systemet kan skrivas på formen d(xi)/dt = inv(M)*(b_tot(t) - A_tot*xi(t)) := g(t, xi), vilket är precis den typ av ekvation ni löste i ALA-B med my_ode! Som ni kanske minns finns det olika metoder för att approximativt lösa system av ODE, t.ex. framåt Euler som vi främst använde i denna kurs. Nu skall vi dock använda metoden bakåt Euler, som numeriskt är mer stabil. Algoritmen för bakåt Euler kan formuleras enligt följande: Inför en partition 0 = t_0 < t_1 < t_2 < ... < t_L = T, med tidssteg dt mellan de diskreta tidpunkterna. Sätt xi^0 = xi(0). För n = 1, 2, 3, ..., L: Beräkna xi^n, som är en approximation till xi(t_n), genom att lösa det linjära ekvationssystemet (M + dt*A_tot)*xi^n = M*xi^(n-1) + dt*b_tot(t_n) Man utgår alltså från begynnelsedata, och beräknar därefter successivt approximationer vid de diskreta tidpunkterna t_1, t_2, ... Dessa approximationer erhålls enligt algoritmen ovan genom att lösa ett linjärt ekvationssystem, vilket alltså måste göras i varje tidssteg. Implementering i Matlab. Viktigt: Se till att samtliga programfiler ligger i samma bibliotek, t.ex. H:, och välj sedan detta bibliotek som Current Directory: Annars finns en risk att Matlab anropar ``fel'' funktioner. Det finns bl.a. en fördefinierad funktion i Matlab som också heter gamma.m, och den vill vi förstås inte skall vara vår permeabilitet... (Om du vill vara säker kan du t.ex. ge kommandot >> which gamma.m så ser du vilken funktion som anropas.) Börja med att skapa en funktionsfil, u0.m, som definierar begynnelsedata: function y = u0(x) % Funktionen u0(x) % % Exempel: y = 0; y = 0; Skapa därefter en funktionsfil, f2.m, som definierar den tidsberoende produktionstermen f(x, t): function source = f2(x, t) % Funktionen f(x, t) % % Exempel: source = (sin(t) + cos(t))*sin(x); source = (sin(t) + cos(t))*sin(x); (Vi kallar filen för f2.m för att inte förväxla den med dess stationära motsvarighet f.m.) Eftersom lastvektorn b = b(t) (potentiellt) beror på tiden kommer denna att behöva räknas om i varje tidssteg. Övriga matriser och vektorer: M, A, Mc, R och rv, är däremot konstanta och kan därmed beräknas en gång för alla innan man går in i tids-loopen. Detta sparar tid och gör programmet effektivare. För att beräkna M kan du använda funktionen MinMassmatrisAssemblerare från Datorstudioövning 2: Kvadratur. L2-projektion. För att beräkna A, Mc, R och rv, kan du i princip använda funktionen MyFirstPoissonAssembler från Datorstudioövning 3: Styvhetsmatris. Robinrandvillkor. Denna funktion räknar dock också ut (den stationära) lastvektorn b vilket inte är önskvärt i detta fall. Kopiera därför MyFirstPoissonAssembler.m till MySecondPoissonAssembler.m, och börja med att ändra funktionshuvudet till function [A, Mc, R, rv] = MySecondPoissonAssembler(p) d.v.s. ta bort b ur argumentlistan (och ändra funktionsnamnet). Radera sedan alla rader som enbart är avsedda för beräkning av b ur koden: ... b = zeros(N, 1); % Ta bort! ... f_xv = f(xv); % Ta bort! f_xh = f(xh); % Ta bort! ... b(i) = b(i) + (f_xv*phi_v_xv + f_xh*phi_v_xh)*h/2; % Ta bort! b(i+1) = b(i+1) + (f_xv*phi_h_xv + f_xh*phi_h_xh)*h/2; % Ta bort! ... Kommentar: Det är i och för sig inte fel att ha dem kvar, men det är ju onödigt och tar extra tid. För att beräkna den tidsberoende lastvektorn b = b(t) modifierar vi istället funktionen MinHLAssemblerare från Datorstudioövning 2: Kvadratur. L2-projektion. Kopiera alltså MinHLAssemblerare.m till t.ex. MinLastAssemblerare.m, och börja med att ändra funktionshuvudet till function b_vek = MinLastAssemblerare(p, tid) d.v.s. lägg till tid till argumentlistan (och ändra funktionsnamnet). Denna funktion måste ha tillgång till aktuell tid, för att kunna skicka denna information vidare till funktionen f2: Ändra anropet av f till ett anrop av f2 enligt ... f_xv = f2(xv, tid); f_xh = f2(xh, tid); ... Kommentar: Det har ingen betydelse att variabeln kallas tid inuti MinLastAssemblerare, men t inuti f2. Det är variabelns värde som skickas och funktionerna kan kalla detta värde för vad de vill. Nu är vi redo att skriva huvudprogrammet. Skapa en fil som du t.ex. kan döpa till MySecondPoissonSolver.m: clear all % Tar bl.a. bort alla variabler från Matlab:s workspace. x_min = 0; % Vänster ändpunkt av rums-intervallet. x_max = pi/2; % Höger ändpunkt av rums-intervallet. p = x_min:(x_max - x_min)/10:x_max; % Definierar rums-meshen. % (10 delintervall.) N = length(p); % Antal nodpunkter. T = pi; % Sluttid. TimeSteps = 20; % Antal tidssteg. tider = 0:T/TimeSteps:T; % Definierar tids-meshen. xi = zeros(N, TimeSteps+1); % Initiera matris där lösningen lagras. % Varje kolonn innehåller lösningens % nodvärden vid en diskret tidpunkt. for i = 1:N % Loopa över alla nodpunkter. x = p(i); % x-koordinat för nod i. xi(i, 1) = u0(x); % Lagra begynnelsenodvärden i kolonn 1. end % Nu assemblerar vi matriser och vektorer som inte beror på tiden: M = MinMassmatrisAssemblerare(p); % Beräkna massmatrisen. [A, Mc, R, rv] = MySecondPoissonAssembler(p); % Beräkna styvhets- % matris, et.c. % Nu startar vi tidsloopen: for n = 1:TimeSteps % Loopa över samtliga tidssteg. n % Skriver ut aktuellt tidssteg. % Kan vara bra att veta hur % långt programmet hunnit! t_h = tider(n+1); % Höger ändpunkt av tidsintervallet. t_v = tider(n); % Vänster ändpunkt av tidsintervallet. dt = t_h - t_v; % Tidsintervallets längd. b = MinLastAssemblerare(p, t_h); % Assemblera lastvektorn vid % tiden t = t_h. (Bakåt Euler!) A_tot = A + Mc + R; b_tot = b + rv; xi_old = xi(:, n); % Lösning vid tiden t = t_v. xi_new = ???; % Bakåt Euler! % Komplettera själv koden! % Använd \ - operatorn för att lösa % det linjära ekvationssystemet. xi(:, n+1) = xi_new; % Lösning vid tiden t = t_h. end % Transponera: (Så fås samma format som Matlabs % egna ode-lösare returnerar.) xi = xi'; tider = tider'; % Nu plottar vi lösningen: mesh(p, tider, xi); xlabel('x') ylabel('t') title('U(x,t)') Testproblem. Testkör nu programmet på följande problem: du/dt - u'' = sin(x)*(sin(t) + cos(t)), 0 < x < pi/2, 0 < t < pi, u(0, t) = u'(pi/2, t) = 0, 0 < t < pi, u(x, 0) = 0, 0 < x < pi/2. Lös t.ex. med 10 delintervall i rummet och 20 delintervall i tiden. Var mycket noga med att kontrollera att samtliga datafiler, a.m o.s.v., är korrekt definierade! Den exakta lösningen till detta problem är (verifiera detta genom insättning i ekvation och randvillkor!) u(x, t) = sin(x)*sin(t) så vi kan testa vår lösning, t.ex. vid tiden t = pi/2: >> MySecondPoissonSolver >> figure % Öppna nytt fönster. >> plot(p, xi(11, :)) % Plottar den approximativa lösningen % vid tiden t = pi/2. Notera att efter % transponering varje rad i xi motsvarar % en diskret tidpunkt. >> hold on >> fplot('sin', [x_min x_max], 'r') % Plottar den exakta lösningen % vid tiden t = pi/2. Notera % skillnaden mellan den exakta % och den beräknade lösningen. % Minns du från ALA-B hur % bakåt Euler beter sig? % Stämmer det i detta fall? Lycka till!!! Tillbaka till kurshemsidan.

Editor: Georgios Foufas Last modified: 2003-04-02