Chalmers


Datorstudioövning 4: Tidsstegningsmetoder

Dagens studioövning består av två delar. (i)Först studerar vi ett enkelt modellproblem: den skalära, ordinära differentialekvationen (ODE),


         du/dt + a(t)*u(t) = f(t),   t > 0,
  (1)
                      u(0) = u0.

Vi skall lösa några olika typ-problem, alla på formen (1), med tre olika tidsstegningsmetoder:

  1. Explicit (framåt) Euler
  2. Implicit (bakåt) Euler
  3. Crank-Nicolson

Målsättningen är att bekanta oss lite med de tre metoderna och lära känna deras styrkor och svagheter. (ii) I den andra delen av dagens övning består uppgiften i att implementera en variant av din PDE-lösare som använder Crank-Nicolson (istället för bakåt Euler, som vi hittills använt) för att lösa systemet av ordinära differentialekvationer (ODE) som uppkommer när vi diskretiserat vår PDE i rummet. Även här skall vi jämföra metoderna genom att lösa några typ-problem. Du skall förstås också behålla din nuvarande bakåt Euler-lösare, och i fortsättningen har vi därmed möjlighet att välja vilken metod vi vill använda. För vissa typer av problem är den ena bättre och vice versa. Det är faktiskt till och med möjligt att använda sig av bägge metoderna för att lösa ett problem, och därmed samtidigt utnyttja bägges goda sidor, vilket vi avslutar med att tala lite om.

Modellproblemet.

Antag att vi önskar lösa (1) för 0 < t < T, och dela därför in tidsintervallet [0 T] i L stycken delintervall:


  0 = t_{0} < t_{1} < ... < t_{n-1} < t_{n} < ... < t_{L} = T

Vi antar för enkelhets skull att varje delintervall, I_{n} = [t_{n-1} t_{n}], har längden dt. Låt U_{n} beteckna en approximation till u(t_{n}). De tre olika tidsstegningsmetoderna lyder nu som ni sett på föreläsningen:

Givet: U_{0} = u(0), beräkna successivt U_{n}, för n = 1,2,...,L, ur,



                   U_{n} = (1 - dt*a_{n-1})*U_{n-1} + dt*f_{n-1}  (framåt Euler)



    (1 + dt*a_{n})*U_{n} = U_{n-1} + dt*f_{n}                      (bakåt Euler)



(1 + 0.5*dt*a_{n})*U_{n} = (1 - 0.5*dt*a_{n-1})*U_{n-1}

                                                                (Crank-Nicolson)
                         + 0.5*dt*(f_{n-1} + f_{n})             


där vi infört beteckningarna a_{n} := a(t_{n}) samt f_{n} := f(t_{n}). Kom också ihåg att metoderna kan härledas genom att integrera (1) över ett tidsintervall, [t_{n-1} t_{n}], samt använda kvadratur för integralberäkningarna. Baserar vi kvadraturen på integrandens värde i vänster ändpunkt fås framåt Euler, höger ändpunkt ger bakåt Euler, och trapetsregeln ger Crank-Nicolson.

För att lösa (1) kan du utgå från följande Matlab-funktion, tre_metoder.m:


  function [] = tre_metoder(u0, T, dt, exaktfinns)

  % u0: begynnelsevärde;  T: sluttid;  dt: tidsstegets längd
  % exaktfinns = 1 <=> exakt lösning känd
  % exaktfinns = 0 <=> exakt lösning okänd

  tidsvektor = [0];     % vi bygger upp en vektor innehållande
                        % de diskreta tidpunkterna
		   
  U_explicit_E = [u0];  % vektor som skall innehålla lösningen
                        % med "framåt Euler"
		     
  U_implicit_E = [u0];  % vektor som skall innehålla lösningen
                        % med "bakåt Euler"     

  U_CN = [u0];          % vektor som skall innehålla lösningen
                        % med "Crank-Nicolson"
		      
  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}
		      
    % framåt Euler:
  
    U_v = U_explicit_E(n);    % U_v = U_{n-1}
    U_h = ...;                % U_h = U_{n};  fyll själv i "metoden"
    U_explicit_E(n+1) = U_h;
  
  
    % bakåt Euler:
  
    U_v = U_implicit_E(n);                     % U_v = U_{n-1}
    U_h = (U_v + dt*f(t_h))/(1 + dt*a(t_h));   % U_h = U_{n}
    U_implicit_E(n+1) = U_h;
  
  
    % Crank-Nicolson:
  
    U_v = U_CN(n);    % U_v = U_{n-1}
    U_h = ...;        % U_h = U_{n};  fyll själv i "metoden"
    U_CN(n+1) = U_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

  % plotta (realdelen av, ifall de är komplexa) lösningarna:

  figure(1)

  plot(tidsvektor, real(U_explicit_E), ':')
  hold on
  plot(tidsvektor, real(U_implicit_E), '--')
  plot(tidsvektor, real(U_CN), '-.')

  if (exaktfinns)  % plotta också den exakta lösningen, om känd
    u_exakt = u(tidsvektor);
    plot(tidsvektor, real(u_exakt), 'g')
  end

  xlabel('t')
  legend('Explicit Euler', 'Implicit Euler', 'Crank-Nicolson', 0)
  hold off


  if (exaktfinns)

  % om den exakta lösningen är känd, plotta (realdelen av) felen:

    figure(2)

    plot(tidsvektor, real(u_exakt - U_explicit_E), ':')
    hold on
    plot(tidsvektor, real(u_exakt - U_implicit_E), '--')
    plot(tidsvektor, real(u_exakt - U_CN), '-.')
    legend('Explicit Euler', 'Implicit Euler', 'Crank-Nicolson', 0)
    title('Fel')
    xlabel('t')
    hold off

  end


  function y = a(t)
  % koefficienten a = a(t)
  y = 1;

  function y = f(t)
  % koefficienten f = f(t)
  y = 0;

  function y = u(t)
  % exakt lösning
  y = exp(-t);

Du får själv fylla i metoderna framåt Euler och Crank-Nicolson. Titta på hur bakåt Euler är gjord.

Lös nu slutligen de tre typ-ekvationerna som gavs på sidorna Tidsstegning 12 - 14, i föreläsningsanteckningarna om Tidsstegning.

Tips: För att lösa ekvationen på sidan Tidsstegning 12 börjar du med att definiera underfunktionerna ("sub functions") a(t) = 1 och f(t) = 0 längst ned i tre_metoder.m (En underfunktion kan enbart anropas av huvudfunktionen, som är samma som namnet på m-filen, och av andra underfunktioner i m-filen.) Eftersom den exakta lösningen är känd i detta fall definierar du den också på samma ställe: u(t) = exp(-t). Anropa nu funktionen så här:


  >> u0 = 1;
  >> T = 1;
  >> dt = 0.1;
  >> exaktfinns = 1;                     % "1" = "sant" i Matlab
  >> tre_metoder(u0, T, dt, exaktfinns)  % Inga utargument finns

Observera att Crank-Nicolson är mera exakt, än de båda "Euler-metoderna".

På sidan Tidsstegning 13 är data komplexa, men det innebär ingen extra svårighet. Matlab räknar med komplexa tal direkt. Definiera endast underfunktionerna a(t) = i och u(t) = cos(t) - i*sin(t). (Verifiera också att detta är en lösning genom att sätta in den i (1).) Sedan fortsätter du som ovan. (Notera att endast realdelar plottas.)

Observera att bakåt Euler verkar "dämpande" medan framåt Euler snarare "förstärker" svängningen. Ånyo är Crank-Nicolson bäst. Nästan för bra, kan man tycka! :) För den exakta lösningen gäller nämligen att |u(t)|^2 = cos(t)^2 + sin(t)^2 = 1, vilket fysikaliskt svarar mot att energin (lägesenergi + kinetisk energi) bevaras. Man kan visa att energin hos lösningen även bevaras för U(t), om vi använder Crank-Nicolson, d.v.s. |U(t)|^2 = konstant = (speciellt) = |U(0)|^2 = |u(0)|^2 = 1, enligt begynnelsevillkoret. Om vi ritar u och U (med t som parameter) i det komplexa talplanet kommer alltså bägge att ligga på enhetscirkeln, och man skulle kunna tro att Crank-Nicolson faktiskt är exakt! Men det är den inte. Zooma in felet i Figure 2 så får du se. Kruxet är att u och U inte befinner sig på samma plats på enhetscirkeln samtidigt: man säger att vi har ett fasfel.

I problemet på sidan Tidsstegning 14 händer något intressant. Crank-Nicolson visar tecken på instabilitet (lösningen oscillerar, men dämpas ändå ut). Framåt Euler blir helt instabil. Endast bakåt Euler ger en god bild av hur den exakta lösningen uppför sig. Det kan vara svårt att se mer än kurvan för framåt Euler ifall man väljer T för stort, eftersom denna lösningskurva efter relativt kort tid blir flera storleksordningar större än de övriga två. Välj då ett litet T, t.ex. 0.2. (Eventuellt ännu mindre men då bör du samtidigt minska steglängden dt.)

Detta är en typisk situation som inträffar när det sker "plötsliga förändringar" i data. Jämför sidan Tidsstegning 19 i föreläsningsanteckningarna. I sådana situationer är därför bakåt Euler att föredra. Å andra sidan är Crank-Nicolson en noggrannare metod och om man kunde skulle man därför helst vilja använda denna. En intressant lösning, t.ex. till problemet på sidan Tidsstegning 14 är följande: ta först ett par steg med bakåt Euler, så att vi "kommer igång" lugnt och fint, byt därefter till Crank-Nicolson! För den som vill kan detta som lite överkurs vara en rolig metod att försöka implementera! Försök att få den stabil, samt mer noggrann än en "ren" bakåt Euler (jämför med den exakta lösningen).

PDE-lösare: Tidsstegning med Crank-Nicolson.

Vi går nu över till reaktions-diffusionsproblemet,


        du/dt - div(a*grad(u)) + c*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 u = u(x,y,t) är den sökta lösningen, och a = a(x,y,t), c = c(x,y,t), f = f(x,y,t), gamma = gamma(x,y,t), g_D = g_D(x,y,t), g_N = g_N(x,y,t), samt u0 = u0(x,y) är data till (2). Antag att vi diskretiserat (2) i rummet och därmed kommit fram till följande system av ODE:

        M*dU/dt + (A + Mc + K)*U = F + G,   t>0,
  (3)
          U(0) = U0,

där U = U(t) är den sökta vektorn (med tidsberoende komponenter/vikter), M är tidsmassmatrisen, A = A(t) är styvhetsmatrisen, Mc = Mc(t) är den "reaktiva" massmatrisen, K = K(t) är randbidragsmatrisen, F = F(t) är lastvektorn, och G = G(t) är randbidragsvektorn. För att härleda Crank-Nicolson för (3) integrerar vi liksom i det skalära fallet över ett tidsintervall, [t_{n-1} t_{n}], samt evaluerar integralerna med trapetsregeln, vilket ger:

   M*(U_{n} - U_{n-1}) +

     0.5*dt*(A_{n-1} + Mc_{n-1} + K_{n-1})*U_{n-1} +

       0.5*dt*(A_{n} + Mc_{n} + K_{n})*U_{n} =

         0.5*dt*(F_{n-1} + G_{n-1} + F_{n} + G_{n})

Vi får därmed följande linjära ekvationssystem som vi måste lösa i varje tidssteg för att beräkna U_{n} givet U_{n-1}:

        ( M + 0.5*dt*(A_{n} + Mc_{n} + K_{n}) )*U_{n} =

  (4)     ( M - 0.5*dt*(A_{n-1} + Mc_{n-1} + K_{n-1}) )*U_{n-1} +

            0.5*dt*(F_{n-1} + G_{n-1} + F_{n} + G_{n})

Övning: Försök nu att implementera en ny version av din lösare som använder (4) istället för bakåt Euler som tidsstegningsmetod. Testa sedan din Crank-Nicolson-lösare på det testproblem som vi använde senast i Datorstudioövning 3: Optimering av kod. Jämför felet med felet du får med bakåt Euler!

Tips: Notera att det inte är så stor skillnad mot bakåt Euler, som du haft förut. I matrisen i vänsterledet är det endast faktorn 0.5 som har tillkommit. I högerledet har det däremot tillkommit ett antal termer som skall evalueras vid tidpunkten t_{n-1}! Vi tar ett exempel: Antag att lastvektorn F beror på tiden men att randvektorn G är konstant. (Vilket t.ex. är fallet i testproblemet från Datorstudioövning 3: Optimering av kod.) Då beräknar du, innan du går in i tidsloopen, G med din "randvektorassemblerare" från förra övningen:


  G = minRandVekAssemblerare(p, e, 0);   % Kan skicka in vad som helst i
                                         % tidsargumentet: används ju ej!

OBS! Namnet på funktionen och argumenten, samt argumentens ordning, behöver förstås ej vara samma i din egen kod.

Inuti tidsloopen måste du sedan anropa din "lastvektorassemblerare" två gånger i varje tidssteg:


  F_v = minLastAssemblerare(p, t, t_v);   % F_v <-> F_{n-1}
  F_h = minLastAssemblerare(p, t, t_h);   % F_h <-> F_{n}

Nu kan du räkna ut den sista raden i (4) genom

  0.5*dt*(F_v + G + F_h + G)  

Vilka matriser som är tidsberoende och vilka som är konstanta varierar förstås från problem till problem, och du får därför vara beredd på att stuva om lite i din kod ibland! Det är ett pris vi får betala för en flexibel och effektiv kod. För ovan nämnda testproblem från förra studioövningen är faktiskt samtliga matriser i vänsterledet konstanta, och du beräknar förstås dessa innan tidsloopen. Notera att du därmed för detta problem med fördel kan använda LU-faktorisering också nu när vi använder Crank-Nicolson!

Lycka till!

Här följer slutligen ett par förslag till ytterligare problem du kan räkna på, men experimentera gärna själv också!

Oscillerande initialdata.

Ett värmeledningsproblem där Omega = [0,1] x [0,1], a = 1, c = f = 0, och u = 0, på randen till Omega. Välj slutligen ett högfrekvent initialdata, t.ex. u0(x,y) = sin(10*pi*x)*sin(10*pi*y). Studera lösningens tidsutveckling. Jämför med bakåt Euler.

Snabb förändring av data.

Samma data som ovan, men välj nu u0(x,y) = 0, f(x,y,t) = 0, 0 < t < 1, t > 2; f(x,y,t) = 25, 1 <= t <= 2. Vad händer vid t=1 och t=2? Jämför med bakåt Euler.

Överkurs: Försök att även i PDE-fallet realisera den kombination av bakåt Euler och Crank-Nicolson som omtalats ovan!

Last modified: Wed Sep 4 12:15:19 MET DST 2002