---

Jeg kan se at emnet har været oppe før her i gruppen, så jeg vil se om ikke
i kan hjælpe mig også.

Jeg sidder også og prøver at få rk4 til at virke med en simulation af
planeters bevægelse omkring hinanden og en central stjerne.

Jeg forstår rk4 som at man beregner k1 som dt * difrentialet af det man
ønsker.
Hvis man ønsker positionen, så er k1 lig dt * farten nu. Det vil i praksis
sige at k1 er ændringen af ens funktionsværdi fra nu til nu+dt.

Imidlertid ser jeg metoden som at man tager et helt tidsskridt med de
nuværende værdier for position,fart og acceleration og ser hvad man ender ud
med. Imidlertid siger formelen dt*den difrentierede funktion af positionen.
Det er mig bekendt farten og kun farten. Der tages altså ikke hensyn til
accelerationen til t=0.

Mit spørgsmål går på om k1 for position skal beregnes som dt*v(t=0) eller
som
dt*v(t=0)+0.5a(t=0)*dt^2 ? Skal accelerationen tages med i beregningen eller
ej?

Andet spørgsmål er om man ikke skal køre rk4 på både position og hastighed,
siden begge afhænger af hinanden.

Skal man først beregne k1 for hastigheden som dt*a(t=0, pos=0) og derefter
for positionen som dt*v(t=0) og i næste skridt sige
k2 for hastigheden er dt*a(t=0.5, pos+0.5k1Pos)

Med andre ord. Skal man beregne k1 for hastigheden og positionen og derefter
bruge v_k1 for at beregne den nye position og pos_k1 for at beregne den nye
acceleration som skal bruges for at beregne v_k2?

Jeg har implementeret det som beskrevet, men jeg kender ikke den analytiske
løsning til et to-legeme problem, så jeg kan kun se at det da ser nogensunde
ud. Jeg ved ikke om det et en helt korrekt implementation.
Hvis en person vil lave en enkelt integration for en planet i kredsløb, så
ville det være en meget stor hjælp. Jeg kunne også godt bruge den analytiske
løsning til problemet, så jeg kan sammenligne forskellige
integrationsmetoder.

På forhånd tak.

---

"Jakob Monstrup" <jms@nomail.thanks> writes:

> Jeg sidder også og prøver at få rk4 til at virke med en simulation
> af planeters bevægelse omkring hinanden og en central stjerne.
[snip]
> Andet spørgsmål er om man ikke skal køre rk4 på både position og
> hastighed, siden begge afhænger af hinanden.

Jeg har ikke selv leget så meget med Runge Kutta, men uden at vide det
vil jeg gætte på at det snarteste er at udregne positionen x(t) og
hastigheden v(t) samtidig. Så slipper du også for problemer ved t=0.

http://www.myphysicslab.com/runge_kutta.html

> Jeg har implementeret det som beskrevet, men jeg kender ikke den analytiske
> løsning til et to-legeme problem, så jeg kan kun se at det da ser nogensunde
> ud.

Ideer:
Prøv at checke om energi og impulsmoment er bevarede.
Prøv evt at gøre den ene masse meget stor og se om du får en elipse.

--
Niels L Ellegaard http://dirac.ruc.dk/~gnalle/

---

> Prøv at checke om energi og impulsmoment er bevarede.
> Prøv evt at gøre den ene masse meget stor og se om du får en elipse.

Det er naturligvis en ide. Energien skal være konstant hvis det kun er to
legemer.
Det kan vel være lige så god en test som sammenligning med den analytiske
løsning.
Det virker velsagtens ikke ved flere legemer, da "slingshot" effekten netop
virker ved at tage energi fra et legeme og overføre til et andet... men der
er så heller ingen analytisk metode der. Takker for ideen, og skammer mig
over ikke selv at have overvejet det

>http://www.myphysicslab.com/runge_kutta.html

Sjovt link og en god forklaring af rk4 for flere variabler. Takker.

---

Jakob Monstrup wrote:
> Jeg sidder også og prøver at få rk4 til at virke med en simulation af
> planeters bevægelse omkring hinanden og en central stjerne.

En simplere løsning er blot at opdatere (i 1D)

v+=a*dt
x+=v*dt

hvor dt løbende bestemmes f.eks. som dt=konst/max(sqrt(a))

Fordi deltax=0.5 a*dt² => dt=konst/sqrt(a)

Du vil så automatisk få højrere præcision når der er store
accelerationer, og derfor brug for højere præcision. Jeg ville
gætte på at denne simple idee vil outperform lavereordens RK
metoder.

> Jeg har implementeret det som beskrevet, men jeg kender ikke den
> analytiske løsning til et to-legeme problem, så jeg kan kun se at det da
> ser nogensunde ud.
Prøv at kigge på den totale energi, hvis den er konstant så er det
godt, hvis ikke så har du problemer.

> Hvis en person vil lave en enkelt integration for en planet i kredsløb, så
> ville det være en meget stor hjælp.

Prøv at se:
http://www.google.com/search?q=integrators+molecular+dynamics&hl=en&lr=&ie=UTF-8&start=10&sa=N


--
Mvh. Carsten Svaneborg
http://www.softwarepatenter.dk

---

In article <gerd2b.gj3.ln@lt30.mpip-mainz.mpg.de>,
zqex@nowhere.on.the.net says...

> En simplere løsning er blot at opdatere (i 1D)
>
> v+=a*dt
> x+=v*dt
>
> hvor dt løbende bestemmes f.eks. som dt=konst/max(sqrt(a))
>
> Fordi deltax=0.5 a*dt² => dt=konst/sqrt(a)
>
> Du vil så automatisk få højrere præcision når der er store
> accelerationer, og derfor brug for højere præcision. Jeg ville
> gætte på at denne simple idee vil outperform lavereordens RK
> metoder.

Forkert - og rigtigt. Det forkerte er at du foreslår euler-integration
(v+=a*dt). Det må man ALDRIG bruge. Det kan bruges som et illustrativt
pædagoisk eksempel, men det er en alt for unøjagtig metode, som i praksis
er ubrugelig.

Det korrekte, du foreslår, er at anvende variabelt tidsstep. Det bør man
altid bruge, når det er muligt, men det er det ikke, hvis man skal
simulere noget i realtid - f.eks. mens resultatet vises som en animation
på skærmen.

Man skal vælge den metode, der kræver færrest mulige regneoperationer for
at opnå den ønskede nøjagtighed. Samtidig skal algoritmen være til at
analysere og implementere. RK4 med variabelt tidsstep betragtes som et
godt kompromis, da den er passende simpel og effektiv.

Med venlig hilsen Sven.

---

Sven Nielsen <sven@SPAMSPAMSPAMscientist.SPAM.com> writes:
> Forkert - og rigtigt. Det forkerte er at du foreslår euler-integration
> (v+=a*dt). Det må man ALDRIG bruge. Det kan bruges som et illustrativt
> pædagoisk eksempel, men det er en alt for unøjagtig metode, som i praksis
> er ubrugelig.

Hmm så vidt jeg kan se var det Verletmetoden han beskrev, for x og v
blev beregnet forskudt, men jeg vil give dig ret i at den nok ikke er
velegnet til simulering af planetbaner. Til gengÃld er metoden meget
anvendt inden for simulering af molekylebevÃgelser. Dette skyldes (sÃ¥
hut jeg visker) at et molekylaert system er meget kaotisk, og saa er
der et eller andet med de fleste metoder bliver omtrendt lige gode,
hvis de opfylder energibevarelse.

Saa vidt jeg husker staar den rigtige forklaring i Hansen og
McDonalds, theory of simple liquids, eller i Allen og Tildesley
computersimulation of simple liquids (eller noget i den retning).

--
Niels L Ellegaard http://dirac.ruc.dk/~gnalle/

---

Hej!

Velocity verlet er symplektisk. Dvs. at udviklingen af
et faserums volumen er en bevaret størrelse.

http://www.ap.univie.ac.at/users/ves/cp0102/dx/node70.html

--
Mvh. Carsten Svaneborg
http://www.softwarepatenter.dk

---

> En simplere løsning er blot at opdatere (i 1D)
>
> v+=a*dt
> x+=v*dt
>
> hvor dt løbende bestemmes f.eks. som dt=konst/max(sqrt(a))
>
> Fordi deltax=0.5 a*dt² => dt=konst/sqrt(a)
>
> Du vil så automatisk få højrere præcision når der er store
> accelerationer, og derfor brug for højere præcision. Jeg ville
> gætte på at denne simple idee vil outperform lavereordens RK
> metoder.

Jeg forstår ikek hvad du mener med 1D er smartere? Smartere som
illustration?
Tanken med at lade tidsskridtet være en funktion af accelerationen er dog
fikst. Har tidligere forsøgt mig med variable tidsskridt, men det har været
ved at tage halve skridt og se på afvigelsen mellem to halve og en hel. Hvis
den var over en given grænse, så blev dt sat ned, ellers op.
Din løsning virker mere simpel og generel.

> Prøv at se:
>
http://www.google.com/search?q=integrators+molecular+dynamics&hl=en&lr=&ie=U
TF-8&start=10&sa=N

Har været der.Flere gange end jeg har lyst til at mindes. Mit problem er at
min matematik er fra gym niveau A + det løse. Jeg kan oftest det jeg skal
bruge, men når det gælder numerisk integration, så har jeg ofte svært ved at
forstå de opgivne formler.
Det vil derfor være en stor hjælp hvis man kan se en step by step løsning af
det aktuelle problem og ikke en videnskablig afhandling.

Imidlertid er der et andet spørgsmål som ikke er besvaret.
Skal jeg beregne ændingen i position som dt*v eller som dt*v + 0.5*a*dt^2 ?
Ifølge formelerne skal jeg for k1 for position beregne dt*position
difrentieret, altså dt*v, men normalt beregnes ænding i position dog med
hensyntagen til a også.

---

Jakob Monstrup wrote:
> Skal man først beregne k1 for hastigheden som dt*a(t=0, pos=0) og derefter
> for positionen som dt*v(t=0) og i næste skridt sige
> k2 for hastigheden er dt*a(t=0.5, pos+0.5k1Pos)
>
> Med andre ord. Skal man beregne k1 for hastigheden og positionen og derefter
> bruge v_k1 for at beregne den nye position og pos_k1 for at beregne den nye
> acceleration som skal bruges for at beregne v_k2?

Det er den rigtige måde at gøre det på. Du skriver ikke andet, så
jeg antager, at du laver en to-dimensionel simulering. For hver
planet skal du derfor til ethvert tidspunkt kende dens position
(x,y) og dens hastighed (vx, vy). Ligningssytemet for én planet er

   x' = vx
   y' = vy
   vx' = ax(x, y, t)
   vy' = ay(x, y, t),

hvor ax(...) og ay(...) er planetens øjeblikkelige acceleration,
som kan beregnes ud fra tyngdeloven og planeternes positioner og
masser. Som det fremgår af to-variabel-eksemplet i linket, Niels
postede, hører der til hver af de dynamiske variabler (x, y, vx,
vy) én "k1-variabel". Disse fire k1'er kunne f.eks. hedde k1x,
k1y, k1vx, k1vy. Som du er inde på i det ovenstående skal først
alle k1'erne beregnes, derefter alle k2'erne ud fra k1'erne, så
alle k3'erne ud fra k2'erne osv (hmm... det lyder som en
tilstandsrapport). Og lad dig i øvrigt ikke forvirre over, at der
helt tilfældigvis er lige så mange sæt k'er, som der er dynamiske
variabler (nemlig fire). Det behøver der ikke at være.

> Jeg har implementeret det som beskrevet, men jeg kender ikke den analytiske
> løsning til et to-legeme problem, så jeg kan kun se at det da ser nogensunde
> ud.

Hvis der er én stjerne og to planeter, er det et
tre-legeme-problem. Det er kaotisk, så løsningen kan generelt ikke
beregnes uden at integrere numerisk.

> Jeg ved ikke om det et en helt korrekt implementation.

Du kan tjekke rigtigheden af rk4-implementationen ved at
sammenligne den med en simpel Euler-integration (hvor x(t+dt) =
x(t) + dt*v(t)). De to skal give (omtrent) samme resultat, når
bare tidsskridtet er så lille, at Euler-algoritmen er rimelig
præcis.

--
Jonas Møller Larsen

---

> Det er den rigtige måde at gøre det på. Du skriver ikke andet, så
> jeg antager, at du laver en to-dimensionel simulering.

i mit testeksempel er det. Reelt er det i 3d, men det kommer jo ud på et.

> Som det fremgår af to-variabel-eksemplet i linket, Niels
> postede, hører der til hver af de dynamiske variabler (x, y, vx,
> vy) én "k1-variabel". Disse fire k1'er kunne f.eks. hedde k1x,
> k1y, k1vx, k1vy. Som du er inde på i det ovenstående skal først
> alle k1'erne beregnes, derefter alle k2'erne ud fra k1'erne, så
> alle k3'erne ud fra k2'erne osv (hmm... det lyder som en
> tilstandsrapport). Og lad dig i øvrigt ikke forvirre over, at der
> helt tilfældigvis er lige så mange sæt k'er, som der er dynamiske
> variabler (nemlig fire). Det behøver der ikke at være.

Mange takker! Det var det jeg gerne ville vide om jeg havde forstået helt
korrekt. Jeg var blevet ret forvirret over vektor notationen y'(t+0.5dt ,
x+0.5*k1) med flere udtryk adskildt af komma. Sædvanligvis er en vektor for
mig tal og ikke underfunktioner.

Der behøver ikke være en k for hver variabel? Du mener hvis en variabel ikke
påvirker beregningen af de andre, så behøver der ikke være en ny k, eller
hvad?
Eller mener du at der ikke behøver være k1,k2,k3,k4(4 stk) for et system med
4 variabler? Hvis det er det sidste, så er jeg dog med så langt Det er
trods alt rk4.

> Hvis der er én stjerne og to planeter, er det et
> tre-legeme-problem. Det er kaotisk, så løsningen kan generelt ikke
> beregnes uden at integrere numerisk.

I mit eksempel nu er det en fast stjerne og en planet. Faktisk flere
planeter, men de påvirker ikke hinanden lige nu. Ville se om jeg kunne få
testet metoden først.

> Du kan tjekke rigtigheden af rk4-implementationen ved at
> sammenligne den med en simpel Euler-integration (hvor x(t+dt) =
> x(t) + dt*v(t)). De to skal give (omtrent) samme resultat, når
> bare tidsskridtet er så lille, at Euler-algoritmen er rimelig
> præcis.

God ide. Særligt når det er to legemer som nu, så er et gennemløb af banen
med lille dt at se som en facitliste. Min rk4 skal så holde sig til samme
bane helt præcis.
Er der dog ikke noget om at rk4 taber energi med tiden?

Der er iøvrigt stadig et uafklaret spørgsmål.
Når jeg skal beregne k1 for hastigheden skal det så være...
k1_v = dt*v + 0.5*a*dt^2 eller
k1_v = dt*v

Som nævnt tidligere, så angiver formelen (som jeg læser den) at det er den
sidste, men samtidigt er det normale jo at man inddrager accelerationen i
sådan en beregning.

---

Jakob Monstrup wrote:
> Mange takker! Det var det jeg gerne ville vide om jeg havde forstået helt
> korrekt. Jeg var blevet ret forvirret over vektor notationen y'(t+0.5dt ,
> x+0.5*k1) med flere udtryk adskildt af komma. Sædvanligvis er en vektor for
> mig tal og ikke underfunktioner.

Der er ikke noget krav om, at input og output til funktioner kun
kan være tal. Inputtet til y' er et tal (tiden) /og/ en vektor med
fire komponenter (planetens postition og hastighed). Så
omformuleret er y' en funktion af i alt fem (tal-)variabler.
Derudover er også y's funktionsværdier 4d vektorer, så i
virkeligheden har vi at gøre med fire selvstændige funktioner af
hver fem variabler. Derfor er vektornotationen rar at have!

> Er der dog ikke noget om at rk4 taber energi med tiden?

RK4 er ikke jo uendelig præcis. Derfor er det ret usandsynligt, at
den bevarer energien eksakt. Men om energien vokser eller aftager,
har jeg ingen idé om. Som andre skriver kan du løbende udregne
energi og impulsmoment som kontrol af præcisionen. Den samlede
impuls skal også være bevaret, men den kan ikke beregnes, så længe
du har en uendelig tung stjerne, som ligger stille (dens impuls er
jo uendelig*0 = ?).

> Der er iøvrigt stadig et uafklaret spørgsmål.
> Når jeg skal beregne k1 for hastigheden skal det så være...
> k1_v = dt*v + 0.5*a*dt^2 eller

Nej. ½*a*dt²-leddet er en del af en analytisk løsning til
problemet, hvor en partikel er udsat for konstant acceleration.
Væk med det.

> k1_v = dt*v

Check.

> Som nævnt tidligere, så angiver formelen (som jeg læser den) at det er den
> sidste,

Brug formelen, dvs højresiden af ligningssystemet i mit forrige
indlæg.

> men samtidigt er det normale jo at man inddrager accelerationen i
> sådan en beregning.

Men ændringen af positionen afhænger jo også (indirekte) af
accelerationen gennem hastigheden, så det er helt, som det skal
være.

--
Jonas Møller Larsen

---

Jonas Møller Larsen wrote:
> Jakob Monstrup wrote:
> > Når jeg skal beregne k1 for hastigheden skal det så være...
> > k1_v = dt*v + 0.5*a*dt^2 eller
[klip]

Hovsa, er du sikker på, at du ikke hellere ville spørge om

> > Når jeg skal beregne k1 for positionen skal det så være...
> > k1_x = dt*v + 0.5*a*dt^2 eller
[klip]

Det er i hvert fald det spørgsmål, jeg svarede på

--
Jonas Møller Larsen

---

> Hovsa, er du sikker på, at du ikke hellere ville spørge om
>
> > > Når jeg skal beregne k1 for positionen skal det så være...
> > > k1_x = dt*v + 0.5*a*dt^2 eller
> [klip]

Ok, så lad gå da Hvis ikke x og v lå så langt fra hinanden på
keyboardet, så ville jeg påstå at det var en slåfejl, men reelt var jeg vist
bare uopmærksom. Du har naturligvis ret.

Tak til dig og i andre for svarene. Nu føler jeg mig lidt mere sikker på
hvad det præcis er der foregår.