ODE (Ordinära DifferentialEkvationer)

Syftet med denna laboration är att Du skall bekanta Dig med några enkla metoder för att lösa ODE-problem. Du kommer också att lära Dig att skriva om ett ODE-problem som ett system av första ordningens ekvationer. Matlabs ode45-rutin används för att beräkna en approximativ lösningen till problemet.

Först en kort repetition från envariabelkursen.
I denna första del kommer Du att beräkna approximativa lösningar till ett begynnelsevärdesproblem:

y' = f(t, y), y(t0) = y0

Problemet kan lösas för hand så att Du enkelt skall kunna beräkna skillnaden mellan approximation och exakt lösning. Detta gör att Du kan  jämföra effektiviteten hos några olika metoder. Metoderna, som Du kommer att testa, är:

  1. Eulers metod; yk+1 = yk + h f(tk, yk)
  2. Heuns metod: yk+1 = yk + 0.5 h [ f(tk, yk) + f(tk+ h, yk+ h f(tk, yk)) ]
  3. En enkel Adams-metod: yk+1 = yk + 0.5 h [ 3 f(tk, yk) - f(tk-1, yk-1) ]

I samtliga fall itererar vi när k = 0, 1, 2, ..., n-1 och vi använder fixt tidssteg h. Så tk = t0 + k h. Begynnelsevärde y0, starttid t0 samt tidssteg h ges nedan. Slut-tiden ges indirekt via h och n. Vi sa under föreläsningen att Eulers metod är alldeles för ineffektiv för att användas i praktiken. Det gäller nu även de andra två metoderna, men de är bättre än Eulers metod, eftersom felet blir mindre (för samma h). I denna lab står "felet" för det globala felet dvs. skillnaden mellan exakt- och approximativt värde. Det finns teori (som vi ej tagit upp) som säger att man kan förvänta sig att felet i Eulers metod är av storleksordningen h, en första ordningens metod. Metod 2 och 3 är andra ordningens metoder och man förväntar sig att felet är av storleksordningen h2. Jag använder ordet "förväntar" eftersom de feluppskattningar som finns är mycket svaga. I denna lab kommer Du att se om förväntningarna infrias. Den metod, ode45, som Du kommer att använda i nästa lab är en kombination av två metoder, en av fjärde- och en av femte ordningen.

Att Heuns metod ger mindre fel, är inte så konstigt. Metoden kräver ju två funktionsberäkningar (beräkningar av funktionen f) i varje iteration. Mer information borde ge en bättre metod. Man programmerar inte Heuns metod som den är skriven ovan. Det finns ju ingen anledning att räkna ut f(tk, yk) två gånger. Heuns metod är en sk enstegsmetod, den utnyttjar bara information som kan räknas fram i det aktuella steget.

Adamsmetoden är en flerstegsmetod, den utnyttjar information från föregående steg (som man ser till att spara). Metoden kräver, liksom Eulers metod, bara en funktionsberäkning i varje steg, men Adamsmetoden ger ett mindre fel. I tillämpningar har man normalt system av ODE, vilket gör att beräkningen av f tar mest datortid. Detta gör att Adamsmetoden blir mycket snabbare än Eulers metod. Man kan skapa Adamsmetoder (liksom andra metoder) av ännu högre ordning, och dessa metoder blir då ännu snabbare.

Ett litet problem med flerstegsmetoder är att de kräver speciell behandling i början av iterationen. När k = 0 har vi ju inget föregående steg (dvs. f(t-1, y-1) saknas). För att lösa det i denna lab, kan Du stega Adamsmetoden för k = 1, 2, ... (dvs. börja med k = 1 och inte k = 0) och utnyttja (t1, y1) från Heuns metod ((t0, y0) känner Du ju).

Använd de tre metoderna på problemet y' = -y + t / 10, t0 = 0, y0 = 1, n = 51, h = 0.1. Beräkna den exakta lösningen och plotta absolutbeloppet av felet, | y(tk) - yk |, k = 1, ..., n med kommandot semilogy. Felen (för de tre metoderna) skall plottas i samma diagram, så att man enkelt kan jämföra. Dra slutsatser! Verkar resonemanget kring metodens ordning och felet i lösningen hålla?



Nu till ett system av ekvationer.
Ett exempel på hur man använder ode45 hittar Du i föreläsningsanteckningarna.

Bakgrund: På en plan yta har vi, på olika ställen, "klistrat fast" en uppsättning laddade partiklar. Vi har dessutom en fri laddad partikel som kan röra sig på planet, och uppgiften går ut på att bestämma den fria partikelns väg (bana). Alla partiklarna har lika stor laddning (och samma tecken).

Lår r(t) vara den fria partikelns position vid tiden t. r(t) är en vektor med två komponenter (en x- och en y-komponent). Vi känner begynnelsepositionen r(0) samt initialhastigheten, r'(0).

Låt p1, ..., pn vara positionerna för de n fixa partiklarna. Under lämpliga villkor (ingen friktion osv) kan rörelse-ekvationen för den fria partikeln skrivas som (jag skriver inte längre ut tiden t):

r'' = c [ (r - p1) / | r - p1 |3 + (r - p2) / | r - p2 |3 + ... + (r - pn) / | r - pn |3 ]

| | betecknar det vanliga avståndet. c är en konstant som har att göra med partikelns massa och kraften mellan partiklarna, använd c = 1e-5. Några fysikaliska detaljer.

Observera att detta är ett andra ordningens system, r är inte en skalär storhet.
 
Skriv om ovanstående problem som ett system av första ordningens ekvationer på formen y' = f(t, y), y(0) = y(0). Din omskrivning skall klara ett godtyckligt värde på n.

Vi vill nu implementera ovanstående omskriving i form av en Matlabfunktion, låt oss kalla den f. Funktionen kommer att anropas (av ode45) med aktuell tidpunkt, tk, och approximation, y(k), och funktionen skall returnera vektorn av derivator, så en lämplig början är något i stil med (tk svarar alltså mot tkoch yk mot y(k)):

function yprim = f(tk, yk)
yprim = [...

tk och yk är ju vanliga variabelnamn, så man kan givetvis skriva sin funktion på följande sätt i stället:

function yprim = f(t, y)
yprim = [...

Matlabfunktionen, f, svarar direkt mot  f i uttrycket y' = f(t, y). Ode-lösaren, ode45, anropar Din funktion för varje tidssteg. Lösaren tillhandahåller (via parametrarna, tk, yk eller t, y) aktuell tid och approximation av lösningen vid denna tidpunkt. Din funktion får använda dessa värden för att beräkna och returnera derivatorna (yprim). Tiden, t, ingår inte explicit i differentialekvationen (däremot beror ju y av t), men Matlab kräver ändå att man har med tiden eftersom ode45 är en allmän rutin där formen på funktionen är bestämd av den som som skrev ode45. Mer om detta.

Använd kolonnvektorer överallt! Det är (som vanligt) ett bekymmer med att skicka extra variabler till funktionen (hjälp med detta). Du kan välja på att anropa ode45 med extra variabler (help ode45), utnyttja global eller att definiera samma variabel på flera ställen.

Tips lagra p1, ..., pn som kolonner i matrisen P. Du kan lätt bestämma n, givet P (t.ex. genom n = size(P, 2);).

ode45 anropas på följande sätt (förutsatt att funktionen heter f och är lagrad på filen f.m).

[t, y] = ode45(@f, [t0:delta_t:t_slut], y0);

t0, delta_t, t_slut samt y0 (kolonnvektorn av begynnelsevärden) måste ha tilldelats värden innan anropet. ode45 returnerar en vektor, t,  av tidpunkter, där vi efterfrågat lösningen, samt en matris y. Varje rad i matrisen svarar mot en tidpunkt och kolonnerna svarar mot de olika lösningskomponenterna i systemet.

t0 är begynnelsetiden (noll i detta exempel) och t_slut är den sista tidpunkten. Eftersom vi senare skall animera rörelseförloppet vill ha utmatning vid bestämda tidpunkter (vi vill inte ha olika tät utmatning, metoden är ju adaptiv, eftersom det kommer att ge felaktiga fartvariationer). Vi styr detta med variabeln delta_t Observera att detta inte är det tidssteg som metoden använder, utan endast steget i utmatningen.
 
Gör ovanstående och testa på fallen nedan. Plotta partikelns bana, inte hastigheten (välj rätt kolonner från y-matrisen med andra ord, y(:, j) är kolonn j. help colon). Använd även mitt lilla animeringsprogram, anim, för att animera rörelseförloppet (hämta denna m-fil och skriv sedan help anim i Matlab). Redovisa koden och rita partikelns färd i x-y-planet. Animeringen går ju inte att lämna in.

  1. Först ett fall med vars hjälp Du kan kontrollera Din funktion.
     
    p1 = (-1, -0.9),   p2 = (1, -0.9),    r(0) = (0.5, -0.9),   r'(0) = (0, 0).
     
    De fixa partiklarna har samma y-koordinat, -0.9, och x-koordinaterna -1 respektive 1. Den fria partikeln får ingen begynnelsefart (noll) och placeras på linjesegmentet mellan de fixa partiklarna. I detta fall skall den fria partikeln svänga, i x-led, med konstant y-koordinat -0.9 (varför?).
     

  2. Detta är ett lite roligare fall, en potentialgrop. Den fria partikeln är inringad av partiklar som ligger som staketstolpar i ett osynligt stängsel. Stängslet utgörs av den växelverkan vi har mellan den fria partikeln och de fixa partiklarna. Den fria partikeln kan komma ur gropen om ger ger den tillräckligt stor begynnelsehastighet. Generera 16 fixa partiklar som ligger på jämna mellanrum på enhetscirkeln, t.ex. genom: 

        e = exp(sqrt(-1) * [0:pi/8:15*pi/8]);
        P = [real(e); imag(e)]; % blir en 2 x 16-matris

    Låt   r(0) = (0.7, 0),   r'(0) = (0, 0.004)

    Testa gärna några andra begynnelsehastigheter, t.ex. (0, 0.0122) respektive (0, 0.009), eller några andra initialplaceringar.

Lämpliga tidsvärden är t0 = 0, delta_t = 10, t_slut = 2000. Om Du har en snabb dator får Du minska delta_t så att partikeln inte rör sig så snabbt vid animationen.

Med mera Matlabprogrammering skulle man givetvis kunna göra ett litet GUI så att man kunde klicka ut positionerna.