Chalmers


Datorstudioövning 6: System av (icke-linjära) ODE

Vi skall idag betrakta följande system av ordinära differentialekvationer,


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

där u(t) = [u_1(t); u_2(t); ...; u_m(t)] och f(t, u) = [f_1(t, u); f_2(t, u); ...; f_m(t, u)] är m x 1-vektorer, och a(t) är en diagonal m x m-matris, där alltså diagonalelementet i rad i motsvarar diffusionskoefficienten för komponent i i PDE-fallet.

komponentform kan (1) skrivas,


        du_1/dt + a_11(t)*u_1(t) = f_1(t, u),   t > 0,
        du_2/dt + a_22(t)*u_2(t) = f_2(t, u),   t > 0,

  (2)   ...

        du_m/dt + a_mm(t)*u_m(t) = f_m(t, u),   t > 0,
  
                          u_i(0) = u0_i.

Notera att ekvationerna kopplas samman genom att f_i, potentiellt, beror av samtliga komponenter i u. (Kommentar: Vi kunde förstås låtit termen a(t)*u(t) ingå i f(t, u), men väljer att skriva ekvationen på formen (1)-(2) för analogi med PDE-fallet.)

Följande funktion, som du kan döpa till fix_ode.m, löser (1) med Crank-Nicolson som tidsstegningsmetod, och använder sig av fixpunktsiteration för att lösa det (potentiellt) icke-linjära ekvationssystemet i varje tidssteg:


  function [tidsvektor, U] = fix_ode(u0, T, dt)

  % Syfte:

  %   Löser ODE-systemet
  %
  %     du/dt + a(t)*u(t) = f(t, u),   t > 0,
  %
  %                  u(0) = u0,
  %
  %   med Crank-Nicolson som tidsstegningsmetod.
  %   Fixpunktsiteration används för lösning av
  %   (icke-linjärt) ekvationssystem i varje tidssteg.
  %   u och f är m x 1 - vektorer. a är en diagonal
  %   m x m - matris.

  % Inargument:

  %   u0: begynnelsevärde; m x 1 - vektor
  %    T: sluttid;
  %   dt: tidsstegets längd

  % Utargument:

  %   tidsvektor: radvektor som innehåller de diskreta tidpunkterna
  %   U:          matris där kolonn n innehåller den approximativa
  %               lösningen (med Crank-Nicolson) vid tiden (n-1)*dt

  % Anropade underfunktioner: (ligger sist i filen)

  %   a:  beräknar "diffusionsmatrisen";   a(t)
  %   f:  beräknar högerledet;   f(t, u)
  
  tidsvektor(1) = 0;  % "begynnelsetid"
  U(:, 1) = u0;       % lägger begynnelsevärde i kolonn 1

  m = size(u0, 1);    % antal ekvationer

  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}
  
    U_v = U(:, n);   % U_v = U_{n-1}; U i vänster ändpunkt av
                     % aktuellt tidsintervall
  
    a_v = feval('a', t_v);
    a_h = feval('a', t_h);		     
    f_v = feval('f', t_v, U_v);
  
  % Fixpunktsiteration:
  
    U_fix = U_v;            % 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 = feval('f', t_h, U_fix);
    
      U_fix_old = U_fix;    % Spar nuvarande "iterat" för att kunna
                            % uppskatta felet när vi beräknat nästa.
    
      for i=1:m             % Loopa över samtliga ekvationer
      
        U_fix(i) = (1/(1 + a_h(i,i)*dt/2))*((1 - a_v(i,i)*dt/2)*U_v(i) + ...
	                                    (f_v(i) + f_h(i))*dt/2);
      
      end
			  
			  
      % Alt.: U_fix = (eye(m) + a_h*dt/2)\((eye(m) - a_v*dt/2)*U_v + ...
      %                                    (f_v + f_h)*dt/2);
    
      fel_fix = norm(U_fix - U_fix_old)    % Uppskattning av fixpunktsfelet.
                                           % Skriv ut.
    
    end
 
    disp('fixpunktsiteration ok')
  
    U_h = U_fix;
  
    U(:, n+1) = U_h;  % Lösning vid tiden t = t_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


  % Underfunktioner:


  function diffusion = a(t)

  diagonal = [1; 2];          % Diagonalelementen i a.

  diffusion = diag(diagonal); % Skapar en diagonal matris med elementen
                              % i "diagonal" på huvuddiagonalen.

  %

  function u_prick = f(t, u)  % u_prick kanske inte är ett så lämpligt
                              % namn, dock, ifall a är skild från 0... :(

  u_prick = zeros(2, 1);      % Initiera högerled.

Studera koden noggrant. Tänk speciellt på några saker:

1. Vi löste redan på förra studiopasset ett system av ordinära differentialekvationer med icke-linjärt högerled, nämligen det som uppkom vid rumsdiskretisering av den bi-stabila ekvationen. Då hade vi alltså en skalär (partiell) differentialekvation och systemet uppstod som sagt vid diskretisering i rummet. På dagens övning tänker vi oss däremot att vi har m stycken (kopplade) ekvationer. (Kommentar: På måndag, då vi skall betrakta system av PDE, får vi kombinera dessa båda "tankesätt"!)

2. Notera hur (under)funktionen f returnerar en vektor där varje komponent motsvarar en ekvation. Också a returnerar (i praktiken) en motsvarande "diffusionsvektor" även om vi valt att lägga den på huvuddiagonalen av en diagonal matris (vilket förenklar notationen i (1)).

3. Betänk mycket noga hur fix-punktsiterationen går till. Skriv själv upp Crank-Nicolsons metod som du får genom att integrera varje ekvation i (2) över ett tidsintervall och approximera integralerna med trapetsregeln. Dividera därefter varje ekvation med (1 + a_ii(t_h)*dt/2) så att de övergår i "fixpunktsformen" u_i = g_i(u). Jämför sedan med koden och observera där hur varje ekvation itereras i tur och ordning. (Kommentar: Det går också att använda Matlab:s matrisnotation och göra allt med ett kommando, se % Alt.: ... om du vill, men att iterera en ekvation i taget påminner mer om hur vi kommer att göra för system av reaktions-diffusionsekvationer i PDE-fallet.)

Testexempel:

Exponentialfunktioner. Testa fallet a_11 = 1, a_22 = 2, f_1 = f_2 = 0. Detta är redan inskrivet i underfunktionerna men kontrollera att du förstår syntaxen! Lös (1) med begynnelsevärde u0 = [1; 1] och plotta lösningen:


  >> u0 = [1; 1];
  >> T = 1;
  >> dt = 0.01;
  >> [tidsvektor, U] = fix_ode(u0, T, dt);
  >> plot(tidsvektor, U)
  >> xlabel('t')
  >> legend('u_1', 'u_2')

Notera att detta är ett linjärt problem varför endast en iteration behövs i varje tidssteg för fixpunktsiterationen att konvergera. Eftersom vi inte heller har någon koppling mellan ekvationerna i detta exempel är det inte svårt att lösa systemet för hand med begynnelsedata u0 = [1; 1] enligt ovan. Gör detta, och jämför sedan den analytiska (exakta) med den beräknade lösningen genom att ge kommandot >> hold on och därefter plotta den exakta lösningen i samma figur (som den beräknade). Tips: Använd tidsvektor som "x-punkter" också för plottning av den analytiska lösningen.

Harar och rävar. Prova att lösa "ODE-motsvarigheten" till populationsproblemet på sidan System 2 i föreläsningsanteckningarna. Byt helt enkelt ut de diffusiva termerna mot +a_i*u_i (med "vår" notation: +a_ii*u_i). Du kan för enkelhets skull sätta a_11 = a_22 = 0, och därefter experimentera med olika "kritiska" parametervärden och initialdata. Gör också en fasplot med kommandot >> plot(U(1,:), U(2,:)). Tips: Inuti f har du tillgång till bägge lösningskomponenterna som u(1) och u(2). Du måste tilldela värden till u_prick(1) och u_prick(2).

"Vågekvationen". (Om du vill.) Prova att skriva om ekvationen u'' + u = 0 som ett system av två första ordningens ekvationer. Lös därefter t.ex. med begynnelsevärdet u(0) = 0; u'(0) = 1. Tips: Jämför med sidorna System 17-19 i föreläsningsanteckningarna. I detta fall blir f = [0; 0], men notera att a inte är diagonal. Du kan därför inte längre iterera en ekvation i taget, utan måste använda alternativet % Alt.: ... som gör "allt på en gång". Jämför med blockmatrisformen på sidan System 20 i föreläsningsanteckningarna.

Last modified: Wed Sep 4 12:16:17 MET DST 2002