|
Datorstudioövning 7: System av PDE
Vi skall idag betrakta system av kopplade
reaktions-diffusionsekvationer på formen:
du_1/dt - div(a_1*grad(u_1)) = f_1,
(1) du_2/dt - div(a_2*grad(u_2)) = f_2, i Omega, t > 0,
...
du_m/dt - div(a_m*grad(u_m)) = f_m,
-a_1*Du_1/Dn = gamma_1*(u_1 - g_D_1) + g_N_1,
(2) -a_2*Du_2/Dn = gamma_2*(u_2 - g_D_2) + g_N_2,
...
-a_m*Du_m/Dn = gamma_m*(u_m - g_D_m) + g_N_m,
på randen till Omega, t > 0,
u_1(x,y,0) = u0_1(x,y),
(3) u_2(x,y,0) = u0_2(x,y), i Omega,
...
u_m(x,y,0) = u0_m(x,y).
Notera att vi måste specificera randvillkor (2) och
begynnelsevillkor (3) till var och en av ekvationerna
i systemet (1), samt att
källtermen i ekvation nummer i,
f_i = f_i(x, y, t, u) = f_i(x, y, t, u_1, u_2, ..., u_m), nu alltså
(potentiellt) beror av samtliga komponenter av lösningen
u = u(x, y, t) = [u_1(x, y, t); u_2(x, y, t); ...; u_m(x, y, t)].
Vi betraktar till att börja med följande exempel med två kopplade ekvationer:
du_1/dt - (2/pi^2)*div(grad(u_1)) = u_1 + u_2,
(4) du_2/dt - (1/pi^2)*div(grad(u_2)) = 0.5*u_2 - u_1,
i Omega = [0, 1] x [0, 1], t > 0,
u_1 = u_2 = 0, 0 < x < 1; y = 1; t > 0,
x = 0; 0 < y < 1; t > 0,
(5)
Du_1/Dn = Du_2/Dn = 0, 0 < x < 1; y = 0; t > 0,
x = 1; 0 < y < 1; t > 0,
(6) u_1(x,y,0) = 0,
u_2(x,y,0) = sin(pi*x/2)*cos(pi*y/2), i Omega.
Den exakta lösningen till (4)-(6) är:
u_1(x,y,t) = sin(t)*sin(pi*x/2)*cos(pi*y/2),
(7)
u_2(x,y,t) = cos(t)*sin(pi*x/2)*cos(pi*y/2).
Verifiera att (7) är en lösning, inklusive randvillkoren (5)
och begynnelsevillkoren (6)!
Kommentar:
Även om vi kommer att ha detta specialfall med två ekvationer som ledstjärna
när vi skriver vårt program, kommer det att fungera precis lika
bra för det allmänna fallet med m stycken ekvationer!
Assemblerare: Våra assemblerare
kommer nu att få lite
mer att stå i! Styvhetsmatrisassembleraren kommer att få assemblera
m stycken styvhetsmatriser (en för varje ekvation),
lastvektorassembleraren kommer att få assemblera m stycken
lastvektorer, och så vidare... Vi går nu igenom våra assemblerare en
i taget och ser vad vi måste förändra:
1. Lastvektorassemblerare. Tanken är nu som sagt att
lastvektorassembleraren skall assemblera m stycken
lastvektorer, en för varje ekvation. (I exemplet (4)-(6) är alltså
m = 2.) För att realisera detta måste vi först
modifiera funktionen f_u_2D.m, som lastvektorassembleraren
anropar vid beräkning av funktionsvärden till källtermen. För systemet (4)
skulle vi t.ex. kunna skriva en funktionsfil, f_sys.m,
enligt:
function source = f_sys(x, y, t, u, delomradesnummer)
% Inargument:
%
% u: 1 x m - matris, där element (1, i) innehåller
% lösningskomponent u_i(x, y, t)
% Utargument:
%
% source: 1 x m - matris, där element (1, i) innehåller
% källtermskomponent f_i(x, y, t, u)
source = zeros(1, 2); % Initiera 1 x 2 - matris.
source(1) = u(1) + u(2); % f_1(x,y,t)
source(2) = 0.5*u(2) - u(1); % f_2(x,y,t)
Tänka noga igenom denna funktion: den returnerar alltså en
radvektor innehållande [f_1 f_2 ... f_m].
Notera också att den förväntar sig att inargumentet
u likaledes är en radvektor innehållande
[u_1 u_2 ... u_m]. Detta måste vi förstås ha i
åtanke när vi skriver vår lastvektorassemblerare.
Vi skriver nu en modifierad version,
minLastAssemblerareUSys.m,
av vår "skalära"
lastvektorassemblerare. Tanken är att den skall returnera
en nnodes x antal_ekv-matris, F, där
kolonn nummer j
innehåller lastvektorn för ekvation nummer j.
(Vi använder fortsättningsvis variabelnamnet antal_ekv
istället för m för att beteckna antalet ekvationer.) I stora
drag skulle den kunna se ut så här:
function F = minLastAssemblerareUSys(p, t, tid, U, antal_ekv);
% U: nnodes x antal_ekv - matris innehållande (nodvärdena av)
% den approximativa lösningens samtliga komponenter vid
% tiden 'tid'. Nodvärdena för komponent u_j skall ligga
% i kolonn nummer j av U.
nnodes = size(p,2); % Antal noder. (= antal kolonner i p)
ntri = size(t,2); % Antal trianglar. (= antal kolonner i t)
F = zeros(nnodes, antal_ekv); % Initiering av lastvektorer.
for n = 1:ntri % loopa över samtliga trianglar
% Ta reda på hörnens nodnummer:
... % som tidigare!
% Räkna ut U:s värden i noderna i, j, och k:
Unodi = U(i,:); % 1 x antal_ekv - matris
Unodj = U(j,:); % 1 x antal_ekv - matris
Unodk = U(k,:); % 1 x antal_ekv - matris
% Räkna ut källtermerna i triangelns hörn vid tiden "tid".
% Notera att fnodi, fnodj, och fnodk blir 1 x antal_ekv - matriser:
fnodi = f_sys(xnodi, ynodi, tid, Unodi, delomradesnummer);
fnodj = f_sys(xnodj, ynodj, tid, Unodj, delomradesnummer);
fnodk = f_sys(xnodk, ynodk, tid, Unodk, delomradesnummer);
% Räkna ut värdena på phi_i, phi_j, och phi_k, i triangelns hörn:
... % som tidigare!
% Hörnkvadratur:
F(i,:) = F(i,:) + ...
(fnodi*phi_inodi + fnodj*phi_inodj + fnodk*phi_inodk)*triangelarea/3;
F(j,:) = F(j,:) + ...
(fnodi*phi_jnodi + fnodj*phi_jnodj + fnodk*phi_jnodk)*triangelarea/3;
F(k,:) = F(k,:) + ...
(fnodi*phi_knodi + fnodj*phi_knodj + fnodk*phi_knodk)*triangelarea/3;
end
Tänk mycket noga igenom hur detta går till. Koden
är som du märker nästan identisk med tidigare, men
eftersom fnodi, fnodj, och fnodk
nu är (rad)vektorer assembleras samtliga ekvationer "på en gång".
(Notera att phi_inodi, ..., och triangelarea
fortfarande är skalärer, men att multiplicera en vektor med en skalär
betyder ju att multiplicera varje komponent i vektorn med skalären.)
Om du är osäker på denna matrissyntax i Matlab så experimentera
lite direkt vid Matlab-prompten, t.ex.:
>> A_test = [1 2; 3 4; 5 6]
>> A_test(2,:) = A_test(2,:) + [2 4]*3
2. Randvektorassemblerare. När du tycker att du förstår
hur lastvektorassembleraren fungerar, kan du på egen hand
försöka modifiera randvektorassembleraren, till t.ex.,
minRandVekAssemblerareSys.m. Du måste förstås också
modifiera randdatafunktionerna så att de
passar till randvillkoren (5). För detta ändamål
behöver du känna till numren på de fyra randsegmenten
för enhetskvadraten på den triangulering vi senare skall räkna på.
Du kan anta att sidan y = 0 har randsegmentnummer 1,
och att de därefter är numrerade 2, 3, 4,
i moturs ordning. Tips: Funktionsfilen som specificerar
permeabiliteten, gamma_sys.m, skulle kunna se ut så här:
function permeabilitet = gamma_sys(x,y,t,rand_nr)
% Funktionen skall returnera en 1 x antal_ekv - matris,
% där element (1, i) innehåller gamma_i.
permeabilitet = zeros(1, 2); % Initiera matris.
if ((rand_nr == 3) | (rand_nr == 4))
permeabilitet(1) = 1e5; % gamma_1(x,y,t)
permeabilitet(2) = 1e5; % gamma_2(x,y,t)
end
Randvektorassembleraren, minRandVekAssemblerareSys.m,
skulle kunna börja så här:
function G = minRandVekAssemblerareSys(p, e, tid, antal_ekv);
nnodes = size(p,2); % Antal noder. (= antal kolonner i p)
nedges = size(e,2); % Antal triangelsidor som ligger på ett
% randsegment. (= antal kolonner i e)
G = zeros(nnodes, antal_ekv); % Initiering av randvektorer.
for l = 1:nedges % Loopa över triangelsidor på
% randen till Omega.
... % osv!
OBS! I denna assemblerare förekommer produkter mellan
gamma och g_D. Eftersom dessa nu är (rad)vektorer
måste du här använda "point-wise"-operationen .*
3. Styvhetsmatrisassemblerare. Börja med att definiera
funktionsfilen a_sys.m:
function diffusion = a_sys(x, y, t, delomradesnummer)
% Funktionen skall returnera en 1 x antal_ekv - matris,
% där element (1, i) innehåller diffusionskoefficienten a_i
diffusion = zeros(1, 2); % Initiera matris.
diffusion(1) = 2/(pi^2); % a_1(x,y,t)
diffusion(2) = 1/(pi^2); % a_2(x,y,t)
När vi nu skall modifiera styvhetsmatrisassembleraren uppkommer
ett problem: vi behöver nu skapa en styvhetsmatris för varje
ekvation. Hur lagrar vi dessa? Jo, vi kan använda oss av Matlab:s
tre-dimensionella matriser! I detta fall skapar vi en matris
A av typen nnodes x nnodes x antal_ekv, där
sedan styvhetsmatrisen för ekvation i fås som
A(:,:,i) som är en "vanlig"
nnodes x nnodes-matris. Experimentera lite:
>> A_test = zeros(4,3,2)
>> A_test(3,2,1) = 5
>> A_test(:,3,2) = [1; 2; 3; 4]
>> A_test(:,:,1)
>> A_test(:,:,2)
>> ...
Skapa nu t.ex. funktionsfilen minStyvhetsAssemblerareSys.m:
function A = minStyvhetsAssemblerareSys(p, t, tid, antal_ekv);
nnodes = size(p,2); % Antal noder. (= antal kolonner i p)
ntri = size(t,2); % Antal trianglar. (= antal kolonner i t)
A = zeros(nnodes, nnodes, antal_ekv); % Initiering av styvhetsmatriser.
for n = 1:ntri % loopa över samtliga trianglar
% Ta reda på hörnens nodnummer:
... % som tidigare!
% Räkna ut diffusionskoefficienterna i triangelns hörn
% vid tiden "tid". Notera att anodi, anodj, och anodk
% är 1 x antal_ekv - matriser.
anodi = a_sys(xnodi, ynodi, tid, delomradesnummer);
anodj = a_sys(xnodj, ynodj, tid, delomradesnummer);
anodk = a_sys(xnodk, ynodk, tid, delomradesnummer);
% Gör om anodi, anodj, och anodk till 1 x 1 x antal_ekv -
% matriser; behövs nedan för att dimensionerna skall stämma
% vid assembleringen. Testa direkt vid Matlab-prompten hur
% detta kommando fungerar, t.ex.,
% >> a_test = [12 9]
% >> a_test = reshape(a_test, 1, 1, 2)
anodi = reshape(anodi, 1, 1, antal_ekv);
anodj = reshape(anodj, 1, 1, antal_ekv);
anodk = reshape(anodk, 1, 1, antal_ekv);
% Räkna ut (de konstanta) gradienterna av phi_i, phi_j,
% och phi_k.
... % som tidigare!
A(i,i,:) = A(i,i,:) + ... % 1 x 1 x antal_ekv - matris!
(anodi*dot(gradphi_i,gradphi_i) + ...
anodj*dot(gradphi_i,gradphi_i) + ...
anodk*dot(gradphi_i,gradphi_i))*triangelarea/3;
... % som tidigare!
end
4. Randmatrisassemblerare. Försök göra denna assemblerare
på egen hand. Börja något i stil med:
function K = minRandMatAssemblerareSys(p, e, tid, antal_ekv);
nnodes = size(p,2); % Antal noder. (= antal kolonner i p)
nedges = size(e,2); % Antal triangelsidor som ligger på ett
% randsegment. (= antal kolonner i e)
K = zeros(nnodes, nnodes, antal_ekv); % Initiering av randmatriser.
for l = 1:nedges % Loopa över triangelsidor på
% randen till Omega.
... % osv.!
Tips: Glöm inte att göra "reshape" av gammanodi
och gammanodj.
5. Tidsmassmatrisassemblerare. Du kan behålla din skalära
assemblerare för tidsmassmatrisen, eftersom den blir densamma för
samtliga ekvationer i (1). (Man kan förstås tänka
sig exempel i vilka man har en koefficient, c_i,
också framför du_i/dt i (1);
i sådana fall måste även denna assemblerare modifieras.)
Lösare:
Slutligen modifierar vi vår lösare, som du t.ex. kan kalla
för MyNLSystemSolver.m:
antal_ekv = 2; % Antal ekvationer.
T = pi/2; % Sluttid.
dt = T/35; % Tidsstegets längd.
k = 0:dt:T;
nnodes = size(p,2);
L = length(k) - 1; % Antal tidssteg.
U = zeros(nnodes, antal_ekv, L+1); % Initiera matris där lösningen lagras.
for i = 1:nnodes
x = p(1,i);
y = p(2,i);
U(i, :, 1) = u0_sys(x,y); % Lagra begynnelsevärden.
end
% Tidsoberoende matriser och vektorer:
M = minTidsMassAssemblerare(p, t);
A_stor = minStyvhetsAssemblerareSys(p, t, 0, antal_ekv);
K_stor = minRandMatAssemblerareSys(p, e, 0, antal_ekv);
G_stor = minRandVekAssemblerareSys(p, e, 0, antal_ekv);
% Beräkna LU-faktoriseringar:
Lower_stor = zeros(nnodes, nnodes, antal_ekv); % Initiera
Upper_stor = zeros(nnodes, nnodes, antal_ekv); % Initiera
for ekv = 1:antal_ekv
A = A_stor(:,:,ekv);
K = K_stor(:,:,ekv);
A_tot = M + 0.5*dt*(A + K);
[Lower, Upper] = lu(A_tot);
Lower_stor(:,:,ekv) = Lower;
Upper_stor(:,:,ekv) = Upper;
end
for n = 1:L % Tidsloop.
tidsintervall = n % Skriver ut aktuellt tidsintervall.
% Kan vara bra att veta hur långt
% programmet hunnit!
t_h = k(n+1); % Höger ändpunkt av tidsintervallet.
t_v = k(n); % Vänster ändpunkt av tidsintervallet.
U_old = U(:, :, n); % Lösning vid tiden t = t_v.
% nnodes x antal_ekv - matris
% Tids- och u-beroende matriser och vektorer:
F_v_stor = minLastAssemblerareUSys(p, t, t_v, U_old, antal_ekv);
% Fixpunktsiteration:
U_fix = U_old; % 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_stor = minLastAssemblerareUSys(p, t, t_h, U_fix, antal_ekv);
U_fix_old = U_fix; % Spar nuvarande "iterat" för att kunna
% uppskatta felet när vi beräknat nästa.
for ekv = 1:antal_ekv
A = A_stor(:,:,ekv); K = K_stor(:,:,ekv);
F_v = F_v_stor(:,ekv); F_h = F_h_stor(:,ekv);
G = G_stor(:,ekv);
Lower = Lower_stor(:,:,ekv); Upper = Upper_stor(:,:,ekv);
hl = (M - 0.5*dt*(A + K))*U_old(:,ekv) + ...
0.5*dt*(F_v + G + F_h + G);
y = Lower\hl;
U_fix(:,ekv) = Upper\y;
end
fel_fix = 0;
for ekv = 1:antal_ekv
fel_fix = max([fel_fix norm(U_fix(:,ekv) - U_fix_old(:,ekv))]);
end
fel_fix % Skriv ut uppskattning av fixpunktsfelet.
end
disp('fixpunktsiteration ok')
U_new = U_fix;
U(:, :, n+1) = U_new; % Lösning vid tiden t = t_h.
end
OBS! Du måste modifiera u0_2D.m till u0_sys.m:
function u0 = u0_sys(x,y)
% Funktionen skall returnera en 1 x antal_ekv - matris, där
% element (1, i) innehåller begynnelsevärde för ekvation 'i'.
u0 = zeros(1, 2); % Initiera matris.
u0(1) = 0; % u0_1(x,y)
u0(2) = sin(pi*x/2)*cos(pi*y/2); % u0_2(x,y)
Testexempel: Prova att lösa
exemplet (4)-(6), t.o.m. sluttiden T = pi/2.
Använd samma mesh som när du räknade på den bi-stabila
ekvationen (hämta den med >> load bistab_a_0_01),
och 35 tidssteg.
Eftersom, se (7),
u_1(x, y, 0) = u_2(x, y, pi/2) och
u_1(x, y, pi/2) = u_2(x, y, 0),
kan du t.ex. jämföra dessa:
>> U1 = reshape(U(:,1,:), nnodes, L+1);
>> U2 = reshape(U(:,2,:), nnodes, L+1);
>> max(abs(U1(:, 1) - U2(:, L+1)))
>> max(abs(U2(:, 1) - U1(:, L+1)))
Du kan förstås också visualisera U1 och U2
var för sig med filmtajm.m, och jämföra med (7).
Filmtajmtips: Sätt lämpligt värde
på tiderperbild (kanske 2 eller 3).
Tänk på att filmtajm.m förväntar sig att "lösningsmatrisen"
som skall visualiseras har ett visst namn: du får ändra på
lämpligt sätt så att det stämmer!
Lycka till!!!
|
