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,

        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,

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!!!

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'

  % ...

  % 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);

  % ...

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);

  % ...

  % ...

  % 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;
     
  % ...

  % ...

  % 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).