---

Jeg spurgte tidligere her i gruppen ang rk4, og fik mange gode svar.
Imidlertid ser min implementation ikke imponerende ud.
Euler med halvt så store tidsskridt som rk4 er faktisk lige så god, så hvis
man ser på beregningstiden så er min rk4 ikke anvendelig.
Jeg appender rk4 delen. Håber en venlig sjæl vil kaste et blik på den og
sige mig hvad det er der er galt. I mine øjne er det korrekt nok, men jeg
tvivler på at jeg har ret i det syn.

rk4_dt er det tidsskridt jeg bruger.
double x,y,vx,vy,ax,ay,m,M,G,maxt; ///sættes i adskildt funktion

for (double t=0;t<maxt;t+=rk4_dt) {
///calculate k1
CalculateAcceleration(x,y);
double k1_x = vx * rk4_dt;
double k1_y = vy * rk4_dt;
double k1_vx = ax * rk4_dt;
double k1_vy = ay * rk4_dt;

///calculate k2
CalculateAcceleration(x+k1_x/2,y+k1_y/2);
double k2_x = (vx+k1_vx*0.5) * rk4_dt;
double k2_y = (vy+k1_vy*0.5) * rk4_dt;
double k2_vx = ax * rk4_dt;
double k2_vy = ay * rk4_dt;

///calculate k3
CalculateAcceleration(x+k2_x/2,y+k2_y/2);
double k3_x = (vx+k2_vx*0.5) * rk4_dt;
double k3_y = (vy+k2_vy*0.5) * rk4_dt;
double k3_vx = ax * rk4_dt;
double k3_vy = ay * rk4_dt;

///calculate k4
CalculateAcceleration(x+k3_x,y+k3_y);
double k4_x = (vx+k3_vx) * rk4_dt;
double k4_y = (vy+k3_vy) * rk4_dt;
double k4_vx = ax * rk4_dt;
double k4_vy = ay * rk4_dt;

x = x + (k1_x + 2 * k2_x + 2 * k3_x + k4_x) / 6;
y = y + (k1_y + 2 * k2_y + 2 * k3_y + k4_y) / 6;
vx = vx + (k1_vx + 2 * k2_vx + 2 * k3_vx + k4_vx) / 6;
vy = vy + (k1_vy + 2 * k2_vy + 2 * k3_vy + k4_vy) / 6;

///tegn x,y
}

void CalculateAcceleration(double x, double y)
{
double dx = 0-x;
double dy = 0-y;
double d = sqrt(dx*dx+dy*dy);
double F = (G*M*m) / (d*d);
double ratio = F/d;
double fx = ratio * dx;
double fy = ratio * dy;

ax = fx/m;
ay = fy/m;
}

Hvis der er noget uafklaret ang koden, så sig endelig til.
Jeg er klar over at det ikke er en programmeringsnewsgroup, men på den anden
side kan jeg ikke forvente at man kender rk4 i en sådan gruppe.

På forhånd tak for enhver hjælp i kan give.

---

Jakob Monstrup wrote:
> Euler med halvt så store tidsskridt som rk4 er faktisk lige så god, så hvis
> man ser på beregningstiden så er min rk4 ikke anvendelig.

Det er ikke nødvendigvis et dårligt tegn, fordi for helt små
tidsskridt vil begge algoritmers baner konvergere mod den eksakte
løsning (dog begrænset af, at computeren kun bruger et endeligt
antal bits (16?) til at repræsentere reelle tal). Måske er det
dét, du ser. Prøv at sætte tidsskridtet op. Dette skal ændre meget
mere på Euler-banen end på RK4-banen, hvis implementationen er
korrekt.

Har du beregningskraft nok til at bruge så små tidsskridt, at
Euler-integrationen er nøjagtig nok(tm), er der ingen grund til at
bruge RK4 (men det har man sædvanligvis ikke).

I praksis kan man vurdere, om simulationen er nøjagtig nok(tm) ved
at se på, om resultaterne er stabile, når tidsskridtet f.eks.
fordobles og halveres.

> Jeg appender rk4 delen. Håber en venlig sjæl vil kaste et blik på den og
> sige mig hvad det er der er galt.

Det ser umiddelbart rigtigt ud, men jeg har kun lige skimmet
koden.

--
Jonas Møller Larsen

---

> Det er ikke nødvendigvis et dårligt tegn, fordi for helt små
> tidsskridt vil begge algoritmers baner konvergere mod den eksakte
> løsning (dog begrænset af, at computeren kun bruger et endeligt
> antal bits (16?) til at repræsentere reelle tal). Måske er det
> dét, du ser. Prøv at sætte tidsskridtet op. Dette skal ændre meget
> mere på Euler-banen end på RK4-banen, hvis implementationen er
> korrekt.

Nah, det er ikke fordi euler er god, at rk4 ikke er bedre.
Hvis jeg vælger et tidsskridt hvor eulers bane er meget ustabil, og giver
rk4 samme skridt, så er den ikke meget bedre end euler.

> Det ser umiddelbart rigtigt ud, men jeg har kun lige skimmet
> koden.

Hm.. ok.. underligt. Synes jo også selv det ser rigtigt ud... bare ikke
resultatet.

---

Jakob Monstrup wrote:
> Nah, det er ikke fordi euler er god, at rk4 ikke er bedre.
> Hvis jeg vælger et tidsskridt hvor eulers bane er meget ustabil, og giver
> rk4 samme skridt, så er den ikke meget bedre end euler.

Det skal den være. Hvis du ikke kan se en markant forskel på
algoritmerne ved et absurd langt tidsskridt på f.eks. en 6.-del af
planetomløbsperioden, må der være en fejl i koden. Angående koden
kan det være et problem, at du har globale og lokale variabler med
samme navn (x og y).

--
Jonas Møller Larsen

---

> Det skal den være. Hvis du ikke kan se en markant forskel på
> algoritmerne ved et absurd langt tidsskridt på f.eks. en 6.-del af
> planetomløbsperioden, må der være en fejl i koden. Angående koden
> kan det være et problem, at du har globale og lokale variabler med
> samme navn (x og y).

Ja, mener også der er fejl. Det er dog ikke de lokale og globale variablers
skyld, da de lokale bruges istedet for de globale som jeg ønsker i mit
eksempel. Måske gør det koden mere uklar, men det er dog ikke der fejlen
putter sig.

Anyway, så takker jeg mange gange for den hjælp du har leveret. Jeg føler
trods alt jeg bedre fårstår rk4 nu. Særligt notationen forvirrede mig
tidligere

---

Jonas Møller Larsen skrev

> Har du beregningskraft nok til at bruge så små tidsskridt, at
> Euler-integrationen er nøjagtig nok(tm), er der ingen grund til at
> bruge RK4 (men det har man sædvanligvis ikke).

Jeg vil gerne lige bemærke, at man sædvanligvis opererer med både
afrundings- og trunkeringsfejl i forbindelse med numerisk integration.
Afrundningsfejl skyldes, at udregningerne sker med endelig præcision
(antal bits), mens trunkeringsfejl skyldes, at man så at sige
trunkerer led over en vis orden i den metode man bruger (Euler er fx.
en første-ordens metode, mens RK4 er en 4-ordens metode).

Kigger man på den samlede fejl for en global løsning til et problem
løst vha numerisk integration med en bestemt metode og
skridt-størrelse, kan man se, at første-ordens metoder stort set er
ubrugelige for selv ganske enkle problemer fordi den lille
skridtstørrelse disse metoder nødvendigvis må bruge giver en meget
stor akkumuleret afrundingsfejl. Hvis man kigger på beregningstiden
for en fast samlet trunkeringsfejl (altså selv hvis man ignorerer
betydningen af afrundningsfejl), så vil højere-ordens metoder ofte
bruge lang mindre tid end lavere-ordens metoder.

Jeg kan på stående fod ikke komme i tanke om et problem, hvor jeg
ville anbefale at man brugte Euler-metoden. Før i tiden ville jeg have
sagt, at Euler måske kunne være god hvis man er tvunget til at have
lav skridtstørrelse og ikke er specielt intereseret i nøjagtighed (fx.
til spil hvor feltet ændrer sig "uforudsigeligt" hele tiden), men jeg
har her i gruppen set hvordan man selv her kan gøre bedre. Desværre er
der erfaringsmæssigt mange der benytter sig af Euler uden at ane hvad
begrænsningerne og konsekvensen af denne metode er.

Som baggrund kan jeg nævne, at vi i kaosgruppen på Fysisk Institut,
DTU, brugte en femte/sjette-ordens indlejret RK metode (RK56) til at
løse langt de fleste ikke-lineære (og ofte kaotiske) problemer. Denne
metode har en lav samlet fejl samtidig med, at den tilbyder
skridtkontrol og interpolering for ganske få ekstra beregninger.


Mvh,
--
Filip Larsen