Reševanje Kepplerjevega enačbe z Runge-Kutta metodo

[apovsic]

Torej imamo enačbo:

a = (GM/r^3)*r = dv/dt

za spremenljivki vzameš v in r in ju potem zapišeš po komponentah(torej x,y,z in njihove ustrezne hitrosti).

Za prvi red je Runge-Kutta(torej f(t,r)):

r(t+h) = r(t) + 6h(k1+2k2+2k3+k4), h je časovni korak.

k1=f(t,r)
k2 = f(t+h/2, r+k1/2)
k3= f(t+h/2, r+k2/2)
k4=f(t+h,r+k3)

Kaj se potem vzame za drugu red? Poskusil sem nadomesti r iz zgornjih enačb z v-ji, vendar se enote niso ujemale. Za r, pa niti ne vem kako bi se lotil.

[apovsic]

Pardon, enote se ujemajo, sem šel še enkrat preverit.

Potem za v uporabim zgornjo enačbo?

Za r pa kar f(r,t) = r*h?

[Aniviller]

No, "red" v smislu reda natancnosti je tukaj itak 4. Ti govoris verjetno o stevilu komponent. Tvoj r,k1,k2,k3,k4 imajo vsi po vec komponent (4, ce si v 2D). Sistem enacb se glasi
\(\frac{d}{dt}\{x,y,v_x,v_y\}=f(\{x,y,v_x,v_y\})=\{v_x,v_y,F_x/m,F_y/m\}\)
kjer je F sila (v tvojem primeru 2 komponenti gravitacijske sile). Enote seveda niso enake za vse komponente "vektorja spremenljivk", se pa ujemajo levo-desno (prvi dve komponenti pravita, da je odvod polozaja hitrost, zadnji dve pa, da je odvod hitrosti pospesek).

Ko vstavljas v k1,k2,k3,k4, vstavljas razlicne x,y,vx,vy v desno stran enacbe (najprej zacetne, potem pa tiste, ki so ze nekoliko premaknjeni s prejsnjimi k-ji).

[apovsic]

V bistvu sem kot red mislil red DE, ker komponente so 2x3(zaradi računanja več objektov, referenčna ravnina je seveda ravnina zemlja-sonce).

Problem se pojavi, ker ne znan razširi uporabe Runge-Kutta iz eno na dvo DE. Torej prehod iz f(t,r) -> f(t,r,v).

[Aniviller]

Ja saj to ti govorim. Skoraj vse numericne metode so za diferencialne enacbe prvega reda. Enacbe drugega reda moras pretvorit na prvi red - tvoj par (r,v) smatraj kot 6-komponentni "vektor", za katerega resujes cel sistem (tudi leva stran ima 6 komponent). Kaj je v tem primeru f sem ti pa zgoraj napisal.

[apovsic]

Oprosti, ampak nič mi ne gre v glavo trenutno. Bolj gledam manj mi je jasno.

[Aniviller]

Vsako diferencialno enacbo visjega reda (2. v tem primeru) lahko pretvoris v prvi red, ce vse odvode razen najvisjega stejes kot neodvisne neznanke. Lahko pokazem na 1D primeru. Ce imas enacbo
\(\ddot{x}=F(x,\dot{x})\)
moras pretvorit na sitem dveh enacb:
\(\dot{x}=v\)
\(\dot{v}=F(x,v)\)
r={x,v} lahko smatras kot neznanko, funkcijo, ki prejema cas in dve komponenti, in vraca dve komponenti, pa imas izrazeno kot {v,F} - sistem je postal dvodimenzionalen, ampak je prvega reda, zato lahko uporabis Runge-Kutta algoritem. Vsakic, ko vidis
f(t,r) je torej misljeno, da v paru {v,F} obe komponenti vzames pri danih vrednostih t in r.

Z drugimi besedami - hitrost je po novem koordinata.

[apovsic]

Aja, itak.

Zdaj je jasno.

Hvala

[Popotnik]

V tem primeru še uporabljam NR (ja, bi šel vsaj na GSL, ampak imam program že dolgo časa spisan) in imam neko teževo. Pač zadeva zbezlja in keč je v parametru eps ali pa v hmin (najmanjši korak). Zdaj, če je hmin res majhen, potem v kodi javi napako, sumim, da zaradi tega, ker je zanj še vedno \(x_{next} = x + hmin\), če je pa eps tudi prenizek, pa rutina manjša manjša korak in že na višji inštanci vrže ven, ker sforsira korak manjši od hmin. Pravzaprav nekako bi najprej rad, da bi program še vedno ločil med dvema doubloma, tudi če je razlika zelo majhna, če je pa kaka bolja rešitev za ta problem, pa toliko bolje.

[Popotnik]

Zdaj, če je hmin res majhen, potem v kodi javi napako, sumim, da zaradi tega, ker je zanj še vedno \(x_{next} = x + hmin\),

tu a-ju (a=1) prištejem b (b=10/b^i).

Koda: Izberi vse

2:      1.1
3:      1.01
4:      1.001
5:      1.0001
6:      1.00001
7:      1.000001
8:      1.0000001
9:      1.00000001
10:     1.000000001
11:     1.0000000001
12:     1.00000000001
13:     1.000000000001
14:     1.0000000000001
15:     1.00000000000001
16:     1.000000000000001
17:     1
18:     1
19:     1
20:     1
21:     1
21-> a+b equals a
22:     1
22-> a+b equals a
23:     1
23-> a+b equals a
24:     1
24-> a+b equals a

Hecno, šele ko imam na 19. decimalnem mestu enico, mi javi, da sta števili enaki.

[Aniviller]

Ja korak itak ne sme bit tako majhen, ce rabis tako majhen korak, da dela, potem je itak metoda neuporabna za tvoj problem. Runge Kutta metode se sicer raje z adaptivnim korakom uporablja. Je to sploh tisti NR ki ima double? Ker marsikatera verzija je kar float, kjer je ze na 7. desetiski decimalki problem.

Za double natancnost ti tam okrog 16. decimalke zmanjka. Ce jaz preverim to, dobim pri 16 ze enakost (tudi ce prav a+b==a primerjas). Je pa tukaj vprasanje kako to preverjas - ce to preverja nekje interno v kodi, potem niti ne ves kaj gleda. Tako da lahko s tem dobi vecjo natancnost. Ni nujno da primerja natanko (a+b)==a. Lahko se tudi enakost preverja interno z vecjo natancnostjo (procesor dela z 80-bitno natancnostjo za notranje operacije, da se izogni zaokrozitvenim napakam, double je 64-bitov).

[Popotnik]

Tako je, uporablja double. In ja, uporablja adaptivni korak. Eps prisili, da korak drastično zmanjša in se pol algoritem zlomi.

[Aniviller]

A ti stvar v sonce pade?
Koliko pa das epsilon? Ker pomen tega epsilona si moras dobro pogledat, po moje ga premajhnega dajes.

[Popotnik]

Ja sj po mojem eps mora biti čimmanjši, ker želim dobiti nek spodoben graf. Kam pade planet točno, niti ne vem, ker nisem risal tega grafa. Imam le graf potencialov v času \(t \to \infty\) v odvisnosti od začetnega pogoja; gre za dva potenciala, ker sta dve telesi (še mimoleteča zvezda polega sonca, gre za znano nalogo najbolj priljubljenega predmeta na ljubljanski šoli fizike ). No grafa teh dveh potencialov sta zelo našpičena, bežita gor in dol; sumim, da je to morda tudi zaradi previsokega eps.

[Aniviller]

V vecini primerov, razen kadar planetek res pade v katero izmed sonc, zadeva ne bi smela biti tako strasno obcutljiva. Pri tej nalogi se tu res ne sme zataknit, je vzorcni primer, generacije in generacije ljudi so z NR prisle lepo skozi. Koliksen pa imas ta epsilon sploh? Ta toleranca (max relativna napaka spremembe pri enem koraku) mora biti toliko, da je stvar po nekem stevilu korakov v mejah zeljene resitve. Kaj manj kot \(10^{-6}\) ni smiselno.