ODE

Syftet med denna laboration är att Du skall lära Dig att skriva om ett ode-problem som ett system av första ordningens ekvationer. Du kommer också att använda Matlabs ode45-rutin  för att beräkna och plotta en approximativ lösningen till problemet. Detta labpm innehåller en hel del text (mest ledning och tips), men det är inte mycket att göra.

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). Om Du vill kan Du tänka Dig magneter istället.

Låt 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) (jag kan inte enkelt skriva r-prick i HTML; i fysik skrivs normalt tidsderivator med en eller flera prickar ovanför funktionen).

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

r'' = c [ (r - p₁) / || r - p₁ ||³ + (r - p₂) / || r - p₂ ||³ + ... + (r - p_n) / || r - p_n ||³ ]

|| || betecknar den vanliga tvånormen (avståndet). c är en konstant som har att göra med partikelns massa och kraften mellan partiklarna, använd c = 1e-5.

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, t_k, och approximation, y^(k), och funktionen skall returnera vektorn av derivator, så en lämplig början är något i stil med:

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

Matlabfunktionen, f, svarar direkt mot  f i uttrycket y' = f(t, y),  tk är t och yk är y. k i tk har jag för att betona att funktionen anropas för varje tidssteg, t_k med approximation y^(k). Man kan givetvis välja andra namn på funktionens parametrar, tid resp. derivator eller t och y, till exempel.
 
Tiden, t, ingår inte explicit i differentialekvationen (däremot beror ju y av t), men Matlab kräver ändå att man har med tk eftersom ode45 är en allmän rutin där formen på funktionen är bestämd av den som skrev ode45. Mer om detta.

Använd kolonnvektorer överallt! Det är (som vanligt) ett bekymmer med att skicka extra variabler till funktionen (se diskussionen i samband med minstakvadratlaborationen). Utnyttja global eller definiera samma variabel på flera ställen.

Tips lagra p₁, ..., p_n 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). Hämta genom att trycka ner högra musknappen på länken och välja Save Link As... (det går inte bra att klippa och klistra in).

1. Först ett fall med vars hjälp Du kan kontrollera Din funktion.
    
   p₁ = (-1, -0.9),   p₂ = (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 = 1, t_slut = 2000. Om det går för snabbt, minska delta_t, om det går för långsamt öka delta_t.

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