|
| numerisk integration... igen Fra : Desilva |
Dato : 01-07-02 13:43 |
|
Nu må jeg lægge enhver eventuel stolthed på hylden og indrømme at jeg på
trods af utallige kodeeksempler og formeler ikke kan forstå runge-kutta
metoden.
Den bliver impelmenteret på direkte forskellige måder rundt omkring. Måder
som ikke engang giver samme resultater, på trods af at de alle burde være
rk4.
Jeg har foran mig en formel sigende
k1=(dx)y'(x,y)
k2=(dx)y'(x+dx/2,y+k1/2)
k3=(dx)y'(x+dx/2,y+k2/2)
k4=(dx)y'(x+dx, y+k3)
y(x+dx)=y(x)+(k1+k2*2+k3*2+k4)/6
Hvis vi siger at y(x) er positionen til tiden x og y'(x) så er hastigheden
til tiden x, hvordan skal
y'(x,y) så læses? Hvad laver y inde i parantesen? k1 skal vel ende ud med at
være hastigheden*dt. Det ville for mig give fin mening hvis der stof y'(x),
men y'(x,y)?? Hvad betyder den notation?
I næste skridt står der så y'(x+dx/2,y+k1/2). Hvad pokker gør jeg her?
Jeg vil umidelbart læse det som at y'(x+dx/2,y+k1/2) betyder hastigheden til
tiden x+0.5dt hvis positionen er pos0+k1/2, men selv med denne kreative
tolkning ender jeg ikke rigtig med noget der kan bruges.
Jeg har indtryk af at denne metode er en variant af midpoint eller improoved
euler, hvor man finder hældningen af "kurven" til t=t+dt*0.5 og derefter
tager hele springet til t+dt med denne hældning, men igen... jeg forstår
ikke notationen her.
Er den en venlig sjæl der vil forbarme sig over min mangel på forståelse og
gøre rede for hvordan y'(x+dx/2,y+k1/2) skal tolkes.
På forhånd mange tak!
| |
Jonas Møller Larsen (01-07-2002)
| Kommentar Fra : Jonas Møller Larsen |
Dato : 01-07-02 16:26 |
|
Desilva wrote:
> k1=(dx)y'(x,y)
> k2=(dx)y'(x+dx/2,y+k1/2)
> k3=(dx)y'(x+dx/2,y+k2/2)
> k4=(dx)y'(x+dx, y+k3)
>
> y(x+dx)=y(x)+(k1+k2*2+k3*2+k4)/6
>
> Hvis vi siger at y(x) er positionen til tiden x og y'(x) så er hastigheden
> til tiden x, hvordan skal
> y'(x,y) så læses?
Det betyder bare, at hastigheden både er en funktion af tiden og af
stedet. Hvis y' til dit formål kun afhænger af x, kan du bare bruge
> k1=(dx)y'(x)
> k2=(dx)y'(x+dx/2)
> k3=(dx)y'(x+dx/2)
> k4=(dx)y'(x+dx)
>
> y(x+dx)=y(x)+(k1+k2*2+k3*2+k4)/6
Jeg ved ikke om dette specialtilfælde af rk4 også har et andet navn. Men
det er jo ækvivalent med
y(x+dx)=y(x) + (y'(x) + 4*y'(x+dx/2) + y'(x+dx))/6.
Jeg tror ikke, du skal håbe på nogen dybere, intuitiv forståelse af rk4,
andet end at man åbenbart beregner dy (ændringen i y) på fire
forskellige måder og til sidst tager et tilpas vægtet gennemsnit af
disse. En udledelse af rk4 starter sikkert med Taylorrækken for y(x+dx)
og kræver, at den totale fejl skal være O(dx^4), hvilket altså er
begrundelsen for, at rk4 ser ud, som den gør. Jeg tror, at intuitionen
let går tabt i den "blackbox" af matematiske manipulationer, som leder
til resultatet.
--
Jonas Møller Larsen
| |
Jonas Møller Larsen (01-07-2002)
| Kommentar Fra : Jonas Møller Larsen |
Dato : 01-07-02 16:41 |
|
Jonas Møller Larsen wrote:
> En udledelse af rk4 starter sikkert med Taylorrækken for y(x+dx)
> og kræver, at den totale fejl skal være O(dx^4)
Dette problem har selvfølgelig flere løsninger - men forskellen mellem
forskellige løsninger (resultaterne af forskellige 4.-ordens algoritmer)
er O(dx^4). Men rk4 udmærker sig ved at være rimelig let at implementere
(5 simple formler i stedet for y(x+dx)=dx*y'(x,y)).
--
Jonas Møller Larsen
| |
Desilva (01-07-2002)
| Kommentar Fra : Desilva |
Dato : 01-07-02 17:09 |
|
> Det betyder bare, at hastigheden både er en funktion af tiden og af
> stedet. Hvis y' til dit formål kun afhænger af x, kan du bare bruge
Så i praksis skal jeg ikke finde en enkelt funktion til beregning af
hastighed baseret på tid og position. Jeg er fri til at tolke "y'(x,y)" som
farten til tiden x, når den tidligere position er y. Jeg ved godt at det er
det samme. Det er mere den generelle tolkning der var problemet før.
Så k2=(dx)y'(x+dx/2,y+k1/2) skal bare læses som dx*hastigheden hvis tiden er
to+0.5dt og positionen er den gamle position+ en halv k1?
I et eksempel med en kugle der falder i et tyngdefelt hvis styrke afhænger
af positionen, skal k2 beregnes som
dt*den hastighed objekt ville have hvis vi går et halvt timestep frem i
tiden positionen er pos0+ en halv k1?
Hvis dette ikke er forkert antaget, så mange tusind tak for at forklare det
så selv jeg forstår det
> Jeg tror ikke, du skal håbe på nogen dybere, intuitiv forståelse af rk4,
> andet end at man åbenbart beregner dy (ændringen i y) på fire
> forskellige måder og til sidst tager et tilpas vægtet gennemsnit af
> disse. En udledelse af rk4 starter sikkert med Taylorrækken for y(x+dx)
> og kræver, at den totale fejl skal være O(dx^4), hvilket altså er
> begrundelsen for, at rk4 ser ud, som den gør. Jeg tror, at intuitionen
> let går tabt i den "blackbox" af matematiske manipulationer, som leder
> til resultatet.
Ja, du har nok ret. Det er bare irreterende at bruge en out-of-the-box
formel uden at forstå hvordan den virker. Midpoint metoden er i det mindste
nem at relatere til virkeligheden.
| |
Jonas Møller Larsen (01-07-2002)
| Kommentar Fra : Jonas Møller Larsen |
Dato : 01-07-02 18:14 |
|
Desilva wrote:
> Så k2=(dx)y'(x+dx/2,y+k1/2) skal bare læses som dx*hastigheden hvis tiden er
> to+0.5dt og positionen er den gamle position+ en halv k1?
Netop.
> I et eksempel med en kugle der falder i et tyngdefelt hvis styrke afhænger
> af positionen, skal k2 beregnes som
> dt*den hastighed objekt ville have hvis vi går et halvt timestep frem i
> tiden positionen er pos0+ en halv k1?
Ja. Bortset fra, at rk4 kun løser første ordens differentialligninger.
Hvis du vil simulere Newtons 2. lov, som er 2. orden i tiden, skal du
først omforme bevægelsesligningen til to koblede førsteordens ligninger:
y'' = F(t,y,y') (N2)
=>
y2' = F(t,y1,y2)
y1' = y2
Funktionen y(t) er udskiftet med y1 = y og y2 = y'. Koblede 1.-ordens
ligninger kan implementeres med formlerne fra rk4 ved at sætte
vektorstreger over alle y'er og k'er og opfatte y1 og y2 som
komponenterne i y-vektor.
> Ja, du har nok ret. Det er bare irreterende at bruge en out-of-the-box
> formel uden at forstå hvordan den virker.
Tja, selvom det måske nogle gange kan synes sådan, så er den sorte
matematikkasse jo ikke magisk. Faktisk kommer der aldrig mere ud af den,
end man putter ind i den. Har du først accepteret input, er der ikke
andet at gøre end også at acceptere output.
--
Jonas Møller Larsen
| |
Desilva (02-07-2002)
| Kommentar Fra : Desilva |
Dato : 02-07-02 19:28 |
|
> Ja. Bortset fra, at rk4 kun løser første ordens differentialligninger.
> Hvis du vil simulere Newtons 2. lov, som er 2. orden i tiden,
Anden orden?... fordi position er fart*tid og fart er acceleration*t.. altså
t^2.
>skal du først omforme bevægelsesligningen til to koblede førsteordens
ligninger:
>
> y'' = F(t,y,y') (N2)
> =>
> y2' = F(t,y1,y2)
> y1' = y2
>
> Funktionen y(t) er udskiftet med y1 = y og y2 = y'. Koblede 1.-ordens
> ligninger kan implementeres med formlerne fra rk4 ved at sætte
> vektorstreger over alle y'er og k'er og opfatte y1 og y2 som
> komponenterne i y-vektor.
Jeg er lidt usikker på hvad du mener her.
Hvis jeg lige fører det tilbage til begreberne fra den faldende kugle, så
siger du vel at
s=0.5*a*t^2 + v0*t +s0 ikke kan løses, men det kan
s=s0+v*t og v=a*t ?
Metoden er herefter at løse begge for hvert skridt?
så k1 beregnes for både s og v i første omgang.
k1_s=dt*(vo) og k1_v=dt*(a0)
k2_s=dt*(v0+0.5k1_s) og k2_v=dt*(a0+0.5k1_a)
Er det korrekt forstået?
Jeg vil lige sige til mit forsvar at jeg sædvanligvis ikke er decideret
langsomt opfattende. Jeg har bare efterhånden ses så mange og forskellige
forklaringer på denne metode at jeg rigtigt bliver forvirret.
| |
Jonas Møller Larsen (02-07-2002)
| Kommentar Fra : Jonas Møller Larsen |
Dato : 02-07-02 23:52 |
|
Desilva wrote:
>>Ja. Bortset fra, at rk4 kun løser første ordens differentialligninger.
>>Hvis du vil simulere Newtons 2. lov, som er 2. orden i tiden,
>
>
> Anden orden?
En n'te ordens differentialligning skal pr. definition indeholde den
søgte funktion differentieret n gange. Ikke flere gange, men gerne
færre. y''' + 2y' = 0 er f.eks. en 3.-ordens differentialligning.
>... fordi position er fart*tid og fart er acceleration*t.. altså
> t^2.
Nærmere d²/dt².
>>skal du først omforme bevægelsesligningen til to koblede førsteordens
>
> ligninger:
>
>>y'' = F(t,y,y') (N2)
>>=>
>>y2' = F(t,y1,y2)
>>y1' = y2
>>
>>Funktionen y(t) er udskiftet med y1 = y og y2 = y'. Koblede 1.-ordens
>>ligninger kan implementeres med formlerne fra rk4 ved at sætte
>>vektorstreger over alle y'er og k'er og opfatte y1 og y2 som
>>komponenterne i y-vektor.
>
>
> Jeg er lidt usikker på hvad du mener her.
> Hvis jeg lige fører det tilbage til begreberne fra den faldende kugle, så
> siger du vel at
> s=0.5*a*t^2 + v0*t +s0 ikke kan løses, men det kan
> s=s0+v*t og v=a*t ?
Prøver du at simulere en faldende kugle, hvor accelerationen er en given
funktion af stedet, tror jeg hurtigt du vil opdage, at du er nødt til at
følge både stedet og hastigheden som funktioner af tiden, for at få det
til at køre. Det er også klart, at kuglens bane er fastlagt, kun når
både sted og hastighed kendes til starttidspunktet. Det har ikke noget
med Euler kontra rk4 at gøre (men med at Newtons lov netop er af 2.
orden - kræver 2 begyndelsesbetingelser). Det giver ikke mening at "løse
s=0.5*a*t^2 + v0*t +s0" (Det er jo i sig selv en løsning af
differentialligningen d²s/dt² = a(s) med konstant acceleration).
Omskrivningen ovenfor svarer til
ds/dt = v
dv/dt = a(s)
- intet dybt i det, men den dobbelte tidsafledede er forsvundet. Med
definitionen y = (s, v) og F = (v, a(s)) kan systemet skrives så simpelt
som dy/dt = F(y) - altså en simpel 1.-ordens differentialligning (dog
med vektorer).
> Metoden er herefter at løse begge for hvert skridt?
> så k1 beregnes for både s og v i første omgang.
Ja, for k1 er en vektor med to komponenter.
> k1_s=dt*(vo) og k1_v=dt*(a0)
> k2_s=dt*(v0+0.5k1_s) og k2_v=dt*(a0+0.5k1_a)
Hvis vi antager, at vi er nået til (s0, v0) i simuleringen, skal
ovenstående være
k1_s=dt*(vo) og k1_v=dt*(a(s0))
k2_s=dt*(v0+0.5k1_v) og k2_v=dt*(a(s0+0.5k1_s))
:
> Er det korrekt forstået?
Mon ikke?
--
Jonas Møller Larsen
| |
Desilva (03-07-2002)
| Kommentar Fra : Desilva |
Dato : 03-07-02 11:14 |
|
> Prøver du at simulere en faldende kugle, hvor accelerationen er en given
> funktion af stedet, tror jeg hurtigt du vil opdage, at du er nødt til at
> følge både stedet og hastigheden som funktioner af tiden, for at få det
> til at køre. Det er også klart, at kuglens bane er fastlagt, kun når
> både sted og hastighed kendes til starttidspunktet. Det har ikke noget
> med Euler kontra rk4 at gøre (men med at Newtons lov netop er af 2.
> orden - kræver 2 begyndelsesbetingelser).
Denne faldende kugles acceleration var konstant. Der var den sposition bare
en funktion af hastigheden som var en funktion af accelerationen og tiden.
> > Metoden er herefter at løse begge for hvert skridt?
> > så k1 beregnes for både s og v i første omgang.
>
> Ja, for k1 er en vektor med to komponenter.
Oh. Da du talte om vektorer tidligere troede jeg du mente ex. hastigheden
opløst i sine x,y og z-komponenter.
Det du mener er istedet at hvis man følger både position og hastighed, så
har man en vektor som indeholder disse to?
> k1_s=dt*(vo) og k1_v=dt*(a(s0))
> k2_s=dt*(v0+0.5k1_v) og k2_v=dt*(a(s0+0.5k1_s))
> > Er det korrekt forstået?
> Mon ikke?
Ja det er jo lige det.
Under alle omstændigheder så mange tak for din tid og dine svar
Det kan være jeg vender tilbage hvis jeg stadig ikke kan få det til at makke
ret.
Måske det er lidt nemmere hvis jeg viser et konkret eksempel, og hvad jeg
prøver at gøre.
Tak igen.
| |
Desilva (01-07-2002)
| Kommentar Fra : Desilva |
Dato : 01-07-02 17:09 |
|
> Det betyder bare, at hastigheden både er en funktion af tiden og af
> stedet. Hvis y' til dit formål kun afhænger af x, kan du bare bruge
Så i praksis skal jeg ikke finde en enkelt funktion til beregning af
hastighed baseret på tid og position. Jeg er fri til at tolke "y'(x,y)" som
farten til tiden x, når den tidligere position er y. Jeg ved godt at det er
det samme. Det er mere den generelle tolkning der var problemet før.
Så k2=(dx)y'(x+dx/2,y+k1/2) skal bare læses som dx*hastigheden hvis tiden er
to+0.5dt og positionen er den gamle position+ en halv k1?
I et eksempel med en kugle der falder i et tyngdefelt hvis styrke afhænger
af positionen, skal k2 beregnes som
dt*den hastighed objekt ville have hvis vi går et halvt timestep frem i
tiden positionen er pos0+ en halv k1?
Hvis dette ikke er forkert antaget, så mange tusind tak for at forklare det
så selv jeg forstår det
> Jeg tror ikke, du skal håbe på nogen dybere, intuitiv forståelse af rk4,
> andet end at man åbenbart beregner dy (ændringen i y) på fire
> forskellige måder og til sidst tager et tilpas vægtet gennemsnit af
> disse. En udledelse af rk4 starter sikkert med Taylorrækken for y(x+dx)
> og kræver, at den totale fejl skal være O(dx^4), hvilket altså er
> begrundelsen for, at rk4 ser ud, som den gør. Jeg tror, at intuitionen
> let går tabt i den "blackbox" af matematiske manipulationer, som leder
> til resultatet.
Ja, du har nok ret. Det er bare irreterende at bruge en out-of-the-box
formel uden at forstå hvordan den virker. Midpoint metoden er i det mindste
nem at relatere til virkeligheden.
| |
|
|