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 p1, ..., pn 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 - p1) / || r - p1 ||3 + (r - p2) / || r - p2 ||3 + ... + (r - pn) / || r - pn ||3 ]
|| || 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, tk, 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,
tk 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 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). 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). |
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.