Test di vari metodi Runge Kutta
Test degli algoritmi Runge Kutta Fehlberg con interpolazione di Horn
Runge Kutta di Cash e Karp con interpolazione cubica
Il problema test è il moto di un pianeta attorno al Sole posto nell'origine. L'asse maggiore dell'ellisse è unitario e l'unico parametro arbitrario è l'eccentricità dell'orbita. Il semiasse minore vale b = √( 1 − e · e ). Qui viene riportata la tabella dell'orbita in funzione della anomalia eccentrica η che determina anche il valore del tempo. Nella pratica il tempo viene usato come variabile indipendente e la anomalia eccentrica η viene dedotta dalla relazione detta equazione di Keplero.
| t | x | y | vx | vy |
|---|
| η − e · sin(η) | cos(η) − e | √( 1 − e · e ) · sin(η) | − sin(η)
————————
( 1 − e · cos(η) ) | √( 1 − e · e ) · cos(η)
——————————
( 1 − e · cos(η) ) |
|---|
!
program xrkf45
use fehlberg
use cashkarp
use testpianeta
implicit none
integer,parameter::q=selected_real_kind(15,300)
real(kind=q)::ec,passorkf
real(kind=q),dimension(0:630,2)::posizioni
integer:: kec,passi,giri,unosi
!
write(unit=*,fmt=*)"Test di vari algoritmi di soluzione di sistemi di eq.diff"
write(unit=*,fmt=*)"Runge Kutta Fehlberg RKF45 con interpolante di Horn"
write(unit=*,fmt=*)"Runge Kutta di Cash-Karp con interpolazione cubica"
write(unit=*,fmt=*)
write(unit=*,fmt=*)"Assegna i millesimi di eccentricita' tra 1 e 999"
read(unit=*,fmt=*) kec
kec=max(1,min(kec,999))
ec=kec
ec=ec/1000
write(unit=*,fmt=*)"Assegna il numero di passi in un dodicesimo di giro"
read(unit=*,fmt=*) passi
passi=max(1,passi)
write(unit=*,fmt=*)"Assegna il numero di giri: stampa l'ultimo"
read(unit=*,fmt=*) giri
write(unit=*,fmt=*)"Inizia Fehlberg..."
call faccio_esempio(ec,passi,giri,posizioni)
passorkf=1
passorkf=(2*atan(passorkf))/(3*passi)
call orbita_esempio(kec,posizioni,passorkf)
write(unit=*,fmt=*)"Inizia Cash Karp..."
call program_test(ec,passi,giri,posizioni)
call orbita_esempio(kec,posizioni,passorkf,"c.htm")
write(unit=*,fmt=*)"Amen - un invio"
read(unit=*,fmt=*)
stop
end program xrkf45
!
L'algoritmo di Fehlberg, anche se non nuovissimo, risulta ancora molto usato quando occorre valutare la posizione di un satellite in istanti prefissati e non dipendenti dal metodo usato per l'integrazione numerica.
L'algoritmo di Fehlberg fornisce due soluzioni, una, quella da usare effettivamente, e l'altra, meno precisa, da usare per controllare il grado di precisione ottenuto. Se la soluzione meno precisa si discosta troppo da quella considerata di riferimento è molto, molto opportuno ridurre il passo di avanzamento nel tempo. In letteratura sono riportati vari criteri per la scelta del passo ma forse semplicemente basta applicare una strategia del genere: se il passo costante provoca errori giudicati eccessivi lo si dimezzi fino a soddisfare i vincoli imposti e quindi, progressivamente, finchè il passo è minore di quello base, si incrementi di un tot percento la lunghezza del passo, passo dopo passo. Ovviamente questo potrebbe far scattare nuovamente il vincolo sull'errore e dunque si instaurerebbe un equilibrio dinamico tra la tendenza a ritornare al passo base e quella a ridurre il passo quanto serve per superare le difficoltà numeriche contingenti.
L'algoritmo di Cash e Karp unito all'interpolazione cubica rappresenta il R.K. più nuovo e consigliabile. Anche in questo caso vengono fornite due soluzioni; la meno precisa serve per il controllo del passo di integrazione.