Datorstudioövning 10: Tidsberoende problem (2D).
Dagens datorstudiouppgift är att inkludera
tidsberoende i din FEM-lösare så att du kan lösa problem
på formen:
du/dt - div(a*grad(u)) + c*u = f(x,y,t), i Omega, 0 < t < T,
(1) -a*du/dn = gamma*(u - g_D) + g_N, på Gamma, 0 < t < T,
u(x,y,0) = u0(x,y), i Omega,
med Robinrandvillkor på randen, Gamma, till Omega. Vi antar liksom
vi gjorde i 1D att f = f(x,y,t) kan bero på
tiden, men att övriga koefficienter
a = a(x,y), c = c(x,y) samt
randvillkor
gamma = gamma(x,y), g_D = g_D(x,y),
g_N = g_N(x,y)
inte ändrar sig med tiden.
Börja med att läsa igenom det utdelade materialet om
Tidsberoende problem (2D).
Gör på egen hand
Exercise 1 - 2.
Att bygga ut din stationära lösare
(MyFirst2DPoissonSolver.m)
till
en tidsberoende lösare
(MySecond2DPoissonSolver.m)
fungerar på precis samma sätt som i 1D.
Det är därför bra att först repetera det du gjorde på
Studioövning 5:
Tidsberoende problem.
Till hjälp följer nedan en kopia (liten font) av den tidsberoende
lösaren i 1D
(MySecondPoissonSolver.m)
samt några tips om vad som
du behöver ändra i 2D.
De delar av koden som det inte sägs något om
kan du låta vara oförändrade. Bara det som du behöver ändra
kommenteras nedan. Skapa alltså en fil
som heter MySecond2DPoissonSolver.m, vilken kan vara en
kopia av MySecondPoissonSolver.m, och gör ändringar
enligt nedan:
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.
Dessa tre rader kan du ta bort. Den första skulle ju radera
matriserna p, e, och t från Matlab:s
workspace!!!
Eventuellt kan du byta ut den mot clear functions, som endast
tar bort kompilerade funktioner från minnet. Rad 2-3 definierade
området i 1D. Detta ritar du ju nu i pdetool.
p = x_min:(x_max - x_min)/10:x_max; % Definierar rums-meshen.
% (10 delintervall.)
Denna rad kan du också ta bort. Nodpunkternas koordinater
finns ju nu
redan i matrisen p som du exporterat från pdetool.
N = length(p); % Antal nodpunkter.
Antalet noder i trianguleringen får du nu genom att räkna
antalet kolonner i p:
N = size(p,2); (Alternativt: nnodes = size(p,2);
för att vara konsistenta med notationen!)
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
Du måste skapa en funktion u0_2D.m som beror av
x och y, samt för varje nod
i beräkna både x och y:
x = p(1,i); y = p(2,i); xi(i, 1) = u0_2D(x,y);
% Nu assemblerar vi matriser och vektorer som inte beror på tiden:
M = MinMassmatrisAssemblerare(p); % Beräkna massmatrisen.
Här anropar du istället funktionen
Min2DMassmatrisAssemblerare.m
från
Studioövning 7:
Kvadratur (2D). Massmatris.:
M = Min2DMassmatrisAssemblerare(p, t);
Notera: Om du på denna studioövning
skrivit om Min2DMassmatrisAssemblerare.m till
en egen variant som använder mittpunktskvadratur kan du
förstås anropa denna istället!
[A, Mc, R, rv] = MySecondPoissonAssembler(p); % Beräkna styvhets-
% matris, et.c.
Precis som i 1D skriver du om
MyFirst2DPoissonAssembler.m
till
MySecond2DPoissonAssembler.m
där endast de konstanta matriserna och vektorerna assembleras.
Du tar alltså bort alla rader som endast används
för beräkning av lastvektorn b. Anropet görs sedan med:
[A, Mc, R, rv] = MySecond2DPoissonAssembler(p, e, t);
% 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!)
Här måste du skriva en funktion som du kan kalla
Min2DLastAssemblerare.m och som anropas med:
b = Min2DLastAssemblerare(p, t, t_h);
Denna funktion skall utföra assembleringen av
lastvektorn b, d.v.s. det som du just plockat
bort ur MyFirst2DPoissonAssembler.m, med den skillnaden
att också aktuell tid, t_h, skickas med som
ett argument och
att anropet av f_2D.m byts ut mot ett anrop av
en funktion f2_2D.m som också kan bero av tiden.
Jämför med hur du hanterade tidsberoendet hos f i 1D.
Ett tips är att använda
Min2DMassmatrisAssemblerare.m som mall:
Ändra funktionshuvudet till
function b = Min2DLastAssemblerare(p, t, t_h);
(Notera att vi inte behöver skicka med ``edge-matrisen''
e som argument eftersom
kurvintegral över randen utförs endast vid beräkning av R
och rv.)
Byt sedan ut raden där M initieras mot en motsvarande
initiering av b, och de nio raderna där
M assembleras mot tre rader där b assembleras.
Lägg slutligen till raderna där f:s nodvärden beräknas,
med tidsberoende adderat:
fnodi = f2_2D(xnodi,ynodi,delomradesnummer,th);
och analogt för noderna j och k.
A_tot = A + Mc + R;
b_tot = b + rv;
xi_old = xi(:, n); % Lösning vid tiden t = t_v.
xi_new = (M + dt*A_tot)\(M*xi_old + dt*b_tot); % Bakåt Euler!
xi(:, n+1) = xi_new; % Lösning vid tiden t = t_h.
end
De återstående raderna kan du ta bort eftersom plottningen är
annorlunda i 2D:
% 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)')
Du kan istället t.ex. använda
kommandot pdesurf:
pdesurf(p, t, xi(:, TimeSteps+1))
xlabel('x')
ylabel('y')
title('U(x,y,t), vid tiden t = T')
och plotta lösningen vid (slut)tiden T. Du kan
förstås också plotta den vid andra tider genom
att byta ut TimeSteps+1 mot någon annan siffra.
Som testproblem kan du använda följande:
du/dt - div(grad(u)) = x*(1-x)*y*(1-y)*cos(t) +
2*(x*(1-x) + y*(1-y))*sin(t),
(2)
i Omega = [0,1] x [0,1], för 0 < t < 2*pi,
u(x,y,t) = 0, på Gamma, för 0 < t < 2*pi,
u(x,y,0) = 0, i Omega.
I detta fall skulle alltså f2_2D.m se ut som följer:
function source = f2_2D(x,y,delomradesnummer,tid)
% Funktionen f(x,y,t)
%
% Exempel: source = 1;
source = x*(1-x)*y*(1-y)*cos(tid) + ...
2*(x*(1-x) + y*(1-y))*sin(tid);
Definiera nu också a_2D.m, c_2D.m,
gamma_2D.m, g_D_2D.m, g_N_2D.m,
och u0_2D.m så att (1) överensstämmer med (2).
Rita och triangulera Omega i pdetool och
exportera matriserna p, e och t.
Notera att
T = 2*pi
och du kan t.ex. använda
TimeSteps = 100;
För att starta programmet skriver du förstås:
>> MySecond2DPoissonSolver
Den exakta lösningen är (verifiera!)
u(x,y,t) = x*(1-x)*y*(1-y)*sin(t)
För t = pi/2 är den
exakta lösningen u(x,y,pi/2) = x*(1-x)*y*(1-y), vilket
du kan jämföra med xi(:, 26).
(Notera: tider(1) = 0,
tider(2) = pi/50, tider(3) = 2*pi/50, ..., tider(26) = 25*pi/50 = pi/2)
Du kan därmed beräkna felet vid tiden t =
pi/2:
>> tid = pi/2; % Använd inte namnet t... :)
>> x = p(1,:)'; % Kolonnvektor.
>> y = p(2,:)';
>> u = x.*(1-x).*y.*(1-y).*sin(tid); % Blir också kolonnvektor...
>> fel = u - xi(:, 26); % ... så dimensionerna
% stämmer överens!
>> pdesurf(p, t, fel) % plottar felet
>> max(abs(fel)) % approximerar max-normen av felet
Jämför gärna din approximativa
lösning med den exakta lösningen vid några fler
tidpunkter.
Lycka till!!!
Tillbaka till kurshemsidan.
|