Kb2 - Differentialekvationer och tekniska beräkningar del A, TMA 205 - 2002/2003

Datorstudioövning 11: 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 9: 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.


Editor: Georgios Foufas
Last modified: 2002-10-07