Datorstudioövning 5: Icke-linjärt högerled
Vi skall idag betrakta följande problem,
du/dt - div(a*grad(u)) = u - u^3, i Omega, t>0,
(1) -a*Du/Dn = 0, på randen till Omega, t>0,
u(x,y,0) = u0(x,y), i Omega,
vilket är ett specialfall av
du/dt - div(a*grad(u)) = f, i Omega, t>0,
(2) -a*Du/Dn = gamma*(u - g_D) + g_N, på randen till Omega, t>0,
u(x,y,0) = u0(x,y), i Omega,
där källtermen, f = f(x, y, t, u), nu också kan bero på
lösningen u. Vi tar inte längre med termen c*u
i vänsterledet
eftersom denna motsvarar (det linjära)
specialfallet f = -c*u.
Innan vi sätter igång med att
modifiera vår kod, med målsättningen att
kunna lösa (icke-linjära)
ekvationer på formen (2),
skall vi tala lite om egenskaperna hos lösningar till (1).
Den bi-stabila ekvationen.
Ekvation (1) kallas för den bi-stabila ekvationen,
eftersom det existerar två stabila, stationära
lösningar, nämligen,
u = -1 och u = 1. Att dessa är lösningar
inses
genom direkt
insättning i (1), varvid vi noterar att
f(u) = u - u^3 = u*(1 - u^2) = u*(1 + u)*(1 - u) = 0, för
bägge dessa värden på u.
Vidare så är du/dt = 0 och grad(u) = (0, 0) (nollvektorn),
eftersom u är konstant;
det senare medför även
att det homogena Neumannrandvillkoret i (1) uppfylls.
Hur vet vi att lösningarna är stabila? Antag
t.ex. att u = 1 och låt sedan u perturberas
(störas) på ett sådant sätt att u >~ 1
(lite större än
1). I detta fall blir f(u) < 0 (utnyttja
faktoriseringen av f ovan) vilket ger ett negativt bidrag till
du/dt, d.v.s., u kommer att sträva tillbaka mot
u = 1. Tänk själv igenom fallet u <~ 1, och
visa att även där u strävar mot att komma tillbaka till
u = 1. Gör nu också en motsvarande stabilitetsanalys
för fallet u = -1!
Kommentar: Det existerar förstås en tredje
stationär lösning till
(1), nämligen, u = 0. Denna är dock inte
stabil. Utsätt t.ex. u för en positiv störning, d.v.s.,
låt u >~ 0. I detta fall blir f(u) > 0 vilket
ger ett positivt bidrag till du/dt, d.v.s.,
u kommer att sträva ännu längre bort från
u = 0. Analogt för små negativa störningar.
Lösningar till (1) kommer därmed typiskt att sträva mot
en av de bägge stabila, stationära
lösningarna u = +/- 1.
Vilken beror av initialdata, u0, på
ett mycket intrikat sätt, och
vägen dit kan vara ganska "krokig". Denna komplexitet hos
lösningar till (1) är karakteristisk för icke-linjära
problem.
Vi börjar med att studera två exempel på lösningar i fallet
då Omega = [0, 1] x [0, 1], och
u0 = cos(2*pi*x*x)*cos(2*pi*y*y).
Lösningarna finns sparade i vårt kursbibliotek,
så börja med att skriva
>> matkfkb för att sätta Matlab:s sökväg.
1. I det första fallet är diffusionskoefficienten a = 0.01,
och initialdata som sagt u0 = cos(2*pi*x*x)*cos(2*pi*y*y).
För att hämta
lösningen till Matlab:s work-space,
skriv
>> load bistab_a_0_01.
Skriver du därefter
>> whos så ser du att du nu har tillgång till
de fem matriserna
p, e, t, dt, och U.
(Jo, dt är en 1x1-matris! :)
I detta fall har beräkningen gjorts på en mesh med 697
noder, och med tidssteget dt = 0.05 fram till sluttiden
T = 20, d.v.s., sammanlagt 400 tidssteg. För att
se meshen,
skriv >> pdemesh(p, e, t). För visualisering av
lösningen, som ligger lagrad kolonnvis i 697 x 401-matrisen
U,
kan du skriva följande script som du t.ex. kan
kalla filmtajm2.m:
% Skapar och visar en film av lösningens tidsutveckling.
% Lösningen skall finnas lagrad kolonnvis i matrisen 'U'.
figure % skapar ett nytt grafikfönster
antaltider = size(U, 2); % antal kolonner i 'U'
tiderperbild = 10; % välj hur ofta du vill
% plotta/"spela in" lösningen:
% 1 <-> vid varje tidssteg
% 2 <-> vid vartannat tidssteg, o.s.v.
antalbilder = ceil(antaltider/tiderperbild); % ceil avrundar uppåt
filmmatris = moviein(antalbilder); % initierar en matris som kan
% lagra angivet antal bilder
pekare = 1; % räknare som anger i vilken kolonn i 'filmmatris'
% nästa bild skall sparas
% Nu loopar vi över önskade tider och "spelar in"
% lösningen vid dessa:
for n = 1:tiderperbild:antaltider
pdeplot(p, e, t, 'xydata', U(:, n), 'colormap', 'hot', 'title', num2str((n-1)*dt))
filmmatris(:, pekare) = getframe; % kopierar aktuellt grafikfönster
% till 'filmmatris'
pekare = pekare + 1;
end
movie(filmmatris, 3) % spelar upp filmen tre gånger!!!
Skriv nu >> filmtajm2 för att spela upp lösningens
tidsutveckling som en film. Notera hur kommandot pdeplot,
med argumentet 'xydata',
används för att skapa en två-dimensionell plot
där färgen anger värdet på
U. Ibland kan detta vara åskådligare
än höjden, som
vi tidigare plottat med pdesurf.
(Kommentar:
Det går att göra mycket mer med pdeplot,
som är mer generellt än pdesurf, se >> help
pdeplot.
Faktum är att pdesurf inte utför något annat
än ett anrop till pdeplot med vissa speciella argument:
t.ex. utför kommandot pdesurf(p, t, U(:, n)) i själva
verket kommandot pdeplot(p, [], t, 'xydata', U(:, n), 'xystyle', 'interp',
'zdata', U(:, n), 'zstyle', 'continuous', 'colorbar', 'off').
Experimentera gärna med pdeplot på egen hand!)
Observera följande karakteristika hos lösningen till
(1): Området Omega tenderar, p.g.a. uttrycket för
f(u),
att delas in i "regioner"
där lösningen är antingen +1 eller -1 med "skikt"
emellan. Dessa regioner "slåss" därefter om "herraväldet".
Diffusionen är förstås också
med och "drar" i lösningen, speciellt i "skikten"
där lösningens krökning är stor. (Man säger att
lösningens rörelse är krökningsdriven.)
I detta fall
ser det ut som om "-1" är på väg att "vinna", även om vi
ännu inte nått ett stationärt tillstånd
vid sluttiden T = 20.
2. I det andra fallet är enda skillnaden att diffusionskoefficienten
a = 0.001, samt att vi räknar fram till sluttiden
T = 30, d.v.s., sammanlagt 600 tidssteg
(dt = 0.05 fortfarande, och samma initialdata,
u0 = cos(2*pi*x*x)*cos(2*pi*y*y), som ovan).
Skriv denna gång >> load bistab_a_0_001.
Med >> whos ser du att du ånyo har tillgång till
de fem matriserna
p, e, t, dt, och U.
Den enda av dessa som har ändrats är alltså U, som nu
är en 697 x 601-matris.
Kör åter >> filmtajm2 och jämför med ovan. Notera att
i detta fall, med lägre diffusion, vi initialt får ännu mer markerade
"regioner", vilka dessutom är mer stillastående än tidigare;
detta beror på att det nu tar längre tid innan diffusionen börjar
inverka.
Det krävs att det bildats "tunna skikt", med stor krökning,
innan termen div(a*grad(u)) kan påverka på allvar
och "driva" lösningen i rumsled.
Inte heller i detta fall har
lösningen hunnit bli stationär vid sluttiden T = 30.
Implementation
Vi skall nu försöka modifiera vår lösare så att den klarar ekvationer
av typ (2). Vi kommer att diskutera tidsstegningsmetoden
Crank-Nicolson,
vilken också är den metod som använts för att beräkna
de bägge lösningarna
ovan. (Tänk dock gärna själv igenom fallet med bakåt Euler,
och modifiera din kod för
denna metod på egen hand.) Vi gör modifikationen i tre steg:
i) Funktionsfil för icke-linjär källterm: "f_u_2D.m"
Först måste vi ändra i funktionsfilen som anropas av
"lastvektorassembleraren". Den skulle t.ex.,
för ekvation (1), kunna se ut så här:
function source = f_u_2D(x, y, t, u, delomradesnummer)
% Funktionen f(x,y,t,u)
%
% Exempel: source = u - u^3;
source = u - u^3;
ii) Lastvektorassemblerare: "minLastAssemblerareU.m"
Assembleraren av lastvektorn behöver,
eftersom källtermen f nu beror på u,
också ta emot aktuell lösning
som ett inargument, och
skulle därför t.ex. kunna ha följande "huvud":
function F = minLastAssemblerareU(p, t, tid, U);
% U: nnodes x 1 - vektor innehållande (nodvärdena av)
% den approximativa lösningen, U, vid tiden 'tid'
Det mesta blir därefter sig likt, men du får lägga till
några rader innan du kan anropa f_u_2D.m:
% ...
% Räkna ut U:s värden i noderna i, j, och k.
Unodi = U(i);
Unodj = U(j);
Unodk = U(k);
% Räkna ut källtermen i triangelns hörn vid tiden "tid":
% (Vi skickar med numret på det delområde i vilket triangeln ligger
% som en parameter, för att på ett enkelt sätt kunna specificera olika
% uttryck för f i skilda delområden.)
fnodi = f_u_2D(xnodi, ynodi, tid, Unodi, delomradesnummer);
fnodj = f_u_2D(xnodj, ynodj, tid, Unodj, delomradesnummer);
fnodk = f_u_2D(xnodk, ynodk, tid, Unodk, delomradesnummer);
% ...
(Kommentar: Detta förutsätter förstås att du använder
hörnkvadratur. Om du använder mittpunktskvadratur får
du tänka efter hur du beräknar U i triangelsidornas
mittpunkter. Tips: U är linjär på varje
triangel.)
iii) Lösare: "MyNLTime2DPoissonSolver.m"
Det som är nytt är dels att vi inte har någon "reaktiv massmatris",
dels att vi måste lösa ett icke-linjärt ekvationssystem i
varje tidssteg. För ekvation (1), där det är
möjligt och önskvärt att
använda LU-faktorisering ifall dt är konstant, fås
i sådana fall (du kan givetvis ha lite andra namn på saker och ting
i din egen kod!):
% ...
% Tidsoberoende matriser och vektorer:
M = minTidsMassAssemblerare(p, t); % Tidsmassmatris
A = minStyvhetsAssemblerare(p, t, 0); % Styvhetsmatris
K = minRandMatAssemblerare(p, e, 0); % "Randmatris"
G = minRandVekAssemblerare(p, e, 0); % "randvektor"
% Beräkna LU-faktorisering:
A_tot = M + 0.5*dt*(A + K); % Crank-Nicolson
[Lower, Upper] = lu(A_tot);
% ...
Vi går nu över till lösandet av ekvationssystemet som tidigare
kan ha sett ut så här:
% ...
% Tidsberoende matriser och vektorer:
F_v = minLastAssemblerare(p, t, t_v);
F_h = minLastAssemblerare(p, t, t_h);
hl = (M - 0.5*dt*(A + Mc + K))*U_old + 0.5*dt*(F_v + G + F_h + G);
y = Lower\hl;
U_new = Upper\y;
% ...
För att nu istället använda fixpunktsiteration, kan du
göra ungefär som följer
(jämför med sidorna Ickelinjära ekvationer 22-23 i
föreläsningsanteckningarna):
% ...
% Tids- och u-beroende matriser och vektorer:
F_v = minLastAssemblerareU(p, t, t_v, U_old);
% Fixpunktsiteration:
U_fix = U_old; % Startgissning.
TOL_fix = dt^3; % Tolerans för Crank-Nicolson.
% (bakåt Euler: dt^2)
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 = minLastAssemblerareU(p, t, t_h, U_fix);
hl = (M - 0.5*dt*(A + K))*U_old + 0.5*dt*(F_v + G + F_h + G);
U_fix_old = U_fix; % Spar nuvarande "iterat" för att kunna
% uppskatta felet när vi beräknat nästa.
y = Lower\hl;
U_fix = Upper\y;
fel_fix = norm(U_fix - U_fix_old) % Uppskattning av fixpunktsfelet.
% Skriv ut.
end
disp('fixpunktsiteration ok')
U_new = U_fix;
% ...
Testexempel: Lös (1) med samma mesh, tidssteg, och
data som i något av
de bägge fallen ovan. (Var noga med att definiera korrekta rand- och initialdata!)
Jämför med "referenslösningen". Tips: Efter att du
skrivit t.ex. >> load bistab_a_0_01 skriver du
>> U_ref = U; och beräknar sedan din egen lösning. Eftersom
det kan ta lite tid att räkna behöver du inte räkna ända fram till
T = 20 (eller 30). Du kan ju ändå jämföra: Säg att
du tagit 25 tidssteg, d.v.s., räknat fram till T =
1.25. Felet vid denna tid får du då som
>> fel = U_ref(:, 26) - U(:, 26); (om vi antar att du kallar
din "lösningsmatris" för U). Uppskatta t.ex. max-normen
av felet: >> max(abs(fel)).
Fler exempel: Prova andra initialdata, t.ex.
u0 = 0.1*cos(2*pi*x*x)*cos(2*pi*y*y) (liten amplitud). Studera
lösningens tidsutveckling med filmtajm2.
(Om du inte kallar din "lösningsmatris" för U måste du
förstås ändra till ditt eget variabelnamn
på motsvarande ställen
i filmtajm2.m, alternativt, döpa om matrisen, t.ex. >> U = xi;
ifall din "lösningsmatris" heter xi, innan du kör filmtajm2.)
Variera diffusionskoefficientens (a:s)
storlek och undersök
hur den påverkar lösningen. Tips:
Använd eventuellt en grövre
mesh om du tycker att beräkningen
tar för lång tid. Välj i så fall dt ~
h (för Crank-Nicolson).
|