Datorstudioövning 6: System av (icke-linjära) ODE
Vi skall idag betrakta följande system av ordinära
differentialekvationer,
du/dt + a(t)*u(t) = f(t, u), t > 0,
(1)
u(0) = u0,
där u(t) = [u_1(t); u_2(t); ...; u_m(t)]
och f(t, u) = [f_1(t, u); f_2(t, u); ...; f_m(t, u)]
är m x 1-vektorer, och a(t) är en
diagonal m x m-matris, där alltså diagonalelementet
i rad i motsvarar diffusionskoefficienten för komponent i
i PDE-fallet.
På komponentform kan (1) skrivas,
du_1/dt + a_11(t)*u_1(t) = f_1(t, u), t > 0,
du_2/dt + a_22(t)*u_2(t) = f_2(t, u), t > 0,
(2) ...
du_m/dt + a_mm(t)*u_m(t) = f_m(t, u), t > 0,
u_i(0) = u0_i.
Notera att ekvationerna kopplas samman genom att f_i,
potentiellt, beror av samtliga komponenter i u.
(Kommentar: Vi kunde förstås låtit termen
a(t)*u(t) ingå i f(t, u), men väljer att skriva ekvationen
på formen (1)-(2) för analogi med PDE-fallet.)
Följande funktion, som du kan döpa till fix_ode.m, löser
(1) med Crank-Nicolson som tidsstegningsmetod, och
använder sig av fixpunktsiteration för att lösa det
(potentiellt) icke-linjära ekvationssystemet i varje tidssteg:
function [tidsvektor, U] = fix_ode(u0, T, dt)
% Syfte:
% Löser ODE-systemet
%
% du/dt + a(t)*u(t) = f(t, u), t > 0,
%
% u(0) = u0,
%
% med Crank-Nicolson som tidsstegningsmetod.
% Fixpunktsiteration används för lösning av
% (icke-linjärt) ekvationssystem i varje tidssteg.
% u och f är m x 1 - vektorer. a är en diagonal
% m x m - matris.
% Inargument:
% u0: begynnelsevärde; m x 1 - vektor
% T: sluttid;
% dt: tidsstegets längd
% Utargument:
% tidsvektor: radvektor som innehåller de diskreta tidpunkterna
% U: matris där kolonn n innehåller den approximativa
% lösningen (med Crank-Nicolson) vid tiden (n-1)*dt
% Anropade underfunktioner: (ligger sist i filen)
% a: beräknar "diffusionsmatrisen"; a(t)
% f: beräknar högerledet; f(t, u)
tidsvektor(1) = 0; % "begynnelsetid"
U(:, 1) = u0; % lägger begynnelsevärde i kolonn 1
m = size(u0, 1); % antal ekvationer
n = 1; % aktuellt tidsintervall
t_v = 0; % vänster ändpunkt i aktuellt tidsintervall,
% d.v.s., t_{n-1}
while t_v < T
t_h = n*dt; % höger ändpunkt i aktuellt tidsintervall,
% d.v.s., t_{n}
U_v = U(:, n); % U_v = U_{n-1}; U i vänster ändpunkt av
% aktuellt tidsintervall
a_v = feval('a', t_v);
a_h = feval('a', t_h);
f_v = feval('f', t_v, U_v);
% Fixpunktsiteration:
U_fix = U_v; % Startgissning.
TOL_fix = dt^3; % Tolerans för Crank-Nicolson.
fel_fix = TOL_fix + 1; % Bara så att villkoret i while-loopen
% nedan uppfylls första gången.
while (fel_fix > TOL_fix)
f_h = feval('f', t_h, U_fix);
U_fix_old = U_fix; % Spar nuvarande "iterat" för att kunna
% uppskatta felet när vi beräknat nästa.
for i=1:m % Loopa över samtliga ekvationer
U_fix(i) = (1/(1 + a_h(i,i)*dt/2))*((1 - a_v(i,i)*dt/2)*U_v(i) + ...
(f_v(i) + f_h(i))*dt/2);
end
% Alt.: U_fix = (eye(m) + a_h*dt/2)\((eye(m) - a_v*dt/2)*U_v + ...
% (f_v + f_h)*dt/2);
fel_fix = norm(U_fix - U_fix_old) % Uppskattning av fixpunktsfelet.
% Skriv ut.
end
disp('fixpunktsiteration ok')
U_h = U_fix;
U(:, n+1) = U_h; % Lösning vid tiden t = t_h.
tidsvektor(n+1) = t_h;
t_v = t_h; % Höger ändpunkt i aktuellt tidsintervall
% blir vänster ändpunkt i nästa tidsintervall.
n = n + 1;
end
% Underfunktioner:
function diffusion = a(t)
diagonal = [1; 2]; % Diagonalelementen i a.
diffusion = diag(diagonal); % Skapar en diagonal matris med elementen
% i "diagonal" på huvuddiagonalen.
%
function u_prick = f(t, u) % u_prick kanske inte är ett så lämpligt
% namn, dock, ifall a är skild från 0... :(
u_prick = zeros(2, 1); % Initiera högerled.
Studera koden noggrant. Tänk speciellt på några saker:
1. Vi löste redan på förra studiopasset ett system av
ordinära differentialekvationer med icke-linjärt högerled,
nämligen det som uppkom vid rumsdiskretisering av den bi-stabila
ekvationen. Då hade vi alltså en skalär (partiell)
differentialekvation och
systemet uppstod som sagt vid diskretisering i rummet. På dagens
övning tänker vi oss däremot att vi har m stycken (kopplade)
ekvationer. (Kommentar:
På måndag, då vi skall betrakta system av PDE, får vi kombinera
dessa båda "tankesätt"!)
2. Notera hur (under)funktionen f returnerar
en vektor där varje komponent motsvarar en ekvation.
Också a returnerar (i praktiken) en motsvarande
"diffusionsvektor" även om vi valt att lägga den på huvuddiagonalen
av en diagonal matris (vilket förenklar notationen i (1)).
3. Betänk mycket noga hur fix-punktsiterationen går till.
Skriv själv upp Crank-Nicolsons metod som du får genom
att integrera varje ekvation i (2) över ett tidsintervall
och approximera integralerna med trapetsregeln. Dividera därefter
varje ekvation med (1 + a_ii(t_h)*dt/2)
så att de övergår
i "fixpunktsformen" u_i = g_i(u). Jämför sedan med koden
och observera där hur varje ekvation itereras i tur och ordning.
(Kommentar: Det går också att använda Matlab:s matrisnotation
och göra allt med ett kommando, se % Alt.: ... om du vill,
men att iterera en ekvation i taget påminner mer om hur vi kommer att
göra för system av reaktions-diffusionsekvationer i PDE-fallet.)
Testexempel:
Exponentialfunktioner. Testa fallet a_11 = 1,
a_22 = 2, f_1 = f_2 = 0. Detta är redan inskrivet
i underfunktionerna men kontrollera att du förstår syntaxen! Lös
(1) med begynnelsevärde u0 = [1; 1] och plotta lösningen:
>> u0 = [1; 1];
>> T = 1;
>> dt = 0.01;
>> [tidsvektor, U] = fix_ode(u0, T, dt);
>> plot(tidsvektor, U)
>> xlabel('t')
>> legend('u_1', 'u_2')
Notera att detta är ett linjärt problem varför
endast en iteration behövs i varje tidssteg för
fixpunktsiterationen att konvergera.
Eftersom vi inte heller har någon koppling mellan ekvationerna i
detta exempel är det inte svårt att lösa systemet för hand
med begynnelsedata u0 = [1; 1] enligt ovan. Gör
detta, och jämför sedan den analytiska (exakta)
med den beräknade lösningen genom att ge
kommandot >> hold on och därefter plotta den exakta lösningen
i samma figur (som den beräknade).
Tips:
Använd tidsvektor som "x-punkter" också för
plottning av den analytiska lösningen.
Harar och rävar. Prova att lösa "ODE-motsvarigheten" till
populationsproblemet
på sidan System 2 i föreläsningsanteckningarna. Byt helt
enkelt ut de diffusiva termerna mot +a_i*u_i
(med "vår" notation: +a_ii*u_i).
Du kan för enkelhets skull sätta a_11 = a_22 = 0, och
därefter experimentera med
olika "kritiska" parametervärden och initialdata.
Gör
också en fasplot med kommandot >> plot(U(1,:), U(2,:)).
Tips:
Inuti f har du tillgång till bägge lösningskomponenterna
som u(1) och u(2).
Du måste tilldela värden till
u_prick(1) och u_prick(2).
"Vågekvationen". (Om du vill.) Prova att skriva om
ekvationen u'' + u = 0 som ett system av två
första ordningens ekvationer. Lös därefter t.ex. med
begynnelsevärdet u(0) = 0; u'(0) = 1. Tips:
Jämför med sidorna System 17-19 i
föreläsningsanteckningarna.
I detta fall blir f = [0; 0], men notera att
a
inte är diagonal. Du kan därför inte längre
iterera en ekvation i taget, utan måste använda
alternativet % Alt.: ... som gör "allt på
en gång". Jämför med blockmatrisformen på
sidan System 20 i föreläsningsanteckningarna.
|