La teoriaRunge-Kutta-Fehlberg RKF4(5)
module fehlberg
! Revisione giugno 2008
implicit none
integer,private,parameter::q=selected_real_kind(15,300)
real(kind=q),dimension(4,11),public::k_esempio,y_esempio
real(kind=q),public:: t_esempio,h_esempio,s_esempio
logical,public::horn_esempio
private:: accelerazioni_esempio
public:: rkf,rkfh,vedoparametri,faccio_esempio
contains
!
! RKF45
!
subroutine vedoparametri()
!
! Questa subroutine non ha finalita' operative ma permette solo
! di verificare che alcune equivalenze tra i coefficienti
! dell'algoritmo di interpolazione e quello di calcolo del passo
! siano effettivamente rispettate.
!
real(kind=q),parameter::bv1=16.0_q/135.0_q,bv3=6656.0_q/12825.0_q, &
bv4=28561.0_q/56430.0_q,bv5=-9.0_q/50.0_q,bv6=2.0_q/55.0_q
!
real(kind=q),parameter::bt1=25.0_q/216.0_q,bt3=1408.0_q/2565.0_q, &
bt4=2197.0_q/4104.0_q,bt5=-1.0_q/5.0_q
!
real(kind=q),parameter::b10=1.0_q,b11=301.0_q/120.0_q,&
b12=-269.0_q/108.0_q,b13=311.0_q/360.0_q
real(kind=q),parameter::b31=7168.0_q/1425.0_q,&
b32=-4096.0_q/513.0_q,b33=14848.0_q/4275.0_q
real(kind=q),parameter::b41=-28561.0_q/8360.0_q, &
b42=199927.0_q/22572.0_q,b43=-371293.0_q/75240.0_q
real(kind=q),parameter::b51=57.0_q/50.0_q,b52=-3.0_q,b53=42.0_q/25.0_q
real(kind=q),parameter::b61=-96.0_q/55.0_q,b62=40.0_q/11.0_q, &
b63=-102.0_q/55.0_q
real(kind=q),parameter::b71=3.0_q/2.0_q,b72=-4.0_q,b73=5.0_q/2.0_q
real(kind=q)::b1,b3,b4,b5,b6,b7
real(kind=q),parameter::s=1.0_q
b1=(b10-s*(b11+s*(b12+s*b13)))
b3=(s*(b31+s*(b32+s*b33)))
b4=(s*(b41+s*(b42+s*b43)))
b5=(s*(b51+s*(b52+s*b53)))
b6=(s*(b61+s*(b62+s*b63)))
b7=(s*(b71+s*(b72+s*b73)))
!
write(unit=*,fmt=*)"In vedoparametri: li confronto sperando che coincidano"
write(unit=*,fmt=*)(bv1+bv3+bv4+bv5+bv6)," deve valere 1"
write(unit=*,fmt=*)(bt1+bt3+bt4+bt5)," deve valere 1"
write(unit=*,fmt=*)" "
write(unit=*,fmt="(a,3f15.7)")" ",bv1,b1,bv1-b1
write(unit=*,fmt="(a,3f15.7)")" ",bv3,b3,bv3-b3
write(unit=*,fmt="(a,3f15.7)")" ",bv4,b4,bv4-b4
write(unit=*,fmt="(a,3f15.7)")" ",bv5,b5,bv5-b5
write(unit=*,fmt="(a,3f15.7)")" ",bv6,b6,bv6-b6
write(unit=*,fmt="(a,3f15.7)")" ",0.0_q,b7,-b7
write(unit=*,fmt=*)"La terza colonna deve avere praticamente tutti zero."
return
end subroutine vedoparametri
!
subroutine rkf(k,y,t,h,horn,accelerazioni)
!
! Le matrici y e k debbono avere le stesse dimensioni.
! La soluzione vecchia è y(:,9) mentre quella nuova e' k(:,9)
! t == tempo all'inizio del passo.
! h == ampiezza del passo
! horn == logical. Se vero occorre prepararsi all'interpolazione
! accelerazioni(n,t,k) == subroutine che calcola l'accelerazione
! dove n indica su che colonna di k va messo il risultato, t rappresenta
! l'istante considerato e k(:,:) e' la matrice di cui si usa la colonna 8.
! In k(:,8) all'uscita c'e' la soluzione meno precisa che serve per
! la stima dell'errore fatto per confronto con la soluzione che e' k(:,9)
! In k(:,10) viene trascritta la soluzione vecchia che serve quando
! occorre interpolare.
!
real(kind=q),intent(in)::t,h
real(kind=q),dimension(1:,1:),intent(in)::y
logical,intent(in)::horn
real(kind=q),dimension(1:,1:),intent(out)::k
interface
subroutine accelerazioni(n,t,k)
integer,intent(in)::n
integer,parameter::q=selected_real_kind(15,300)
real(kind=q),intent(in)::t
real(kind=q),intent(in out),dimension(1:,1:)::k
end subroutine accelerazioni
end interface
real(kind=q),parameter::c2=1.0_q/4.0_q,c3=3.0_q/8.0_q, &
c4=12.0_q/13.0_q,c5=1.0_q,c6=1.0_q/2.0_q
real(kind=q),parameter::a21=1.0_q/4.0_q
real(kind=q),parameter::a31=3.0_q/32.0_q,a32=9.0_q/32.0_q
real(kind=q),parameter::a41=1932.0_q/2197.0_q, &
a42=-7200.0_q/2197.0_q,a43=7296.0_q/2197.0_q
real(kind=q),parameter::a51=439.0_q/216.0_q, &
a52=-8.0_q,a53=3680.0_q/513.0_q,a54=-845.0_q/4104.0_q
real(kind=q),parameter::a61=-8.0_q/27.0_q,a62=2.0_q, &
a63=-3544.0_q/2565.0_q,a64=1859.0_q/4104.0_q,a65=-11.0_q/40.0_q
real(kind=q),parameter::bt1=25.0_q/216.0_q,bt3=1408.0_q/2565.0_q, &
bt4=2197.0_q/4104.0_q,bt5=-1.0_q/5.0_q
real(kind=q),parameter::bv1=16.0_q/135.0_q,bv3=6656.0_q/12825.0_q, &
bv4=28561.0_q/56430.0_q,bv5=-9.0_q/50.0_q,bv6=2.0_q/55.0_q
real(kind=q),parameter::h1=1.0_q/6.0_q,h5=1.0_q/6.0_q,h6=2.0_q/3.0_q
!
k(:,8) =y(:,9)
call accelerazioni(1,t,k)
k(:,8) =y(:,9)+h*a21*k(:,1)
call accelerazioni(2,t+c2*h,k)
k(:,8) =y(:,9)+h*(a31*k(:,1)+a32*k(:,2))
call accelerazioni(3,t+c3*h,k)
k(:,8) =y(:,9)+h*(a41*k(:,1)+a42*k(:,2)+a43*k(:,3))
call accelerazioni(4,t+c4*h,k)
k(:,8) =y(:,9)+h*(a51*k(:,1)+a52*k(:,2)+a53*k(:,3)+a54*k(:,4))
call accelerazioni(5,t+c5*h,k)
k(:,8) =y(:,9)+h*(a61*k(:,1)+a62*k(:,2)+a63*k(:,3)+a64*k(:,4)+a65*k(:,5))
call accelerazioni(6,t+c6*h,k)
if(horn) then
k(:,8)=y(:,9)+h*(h1*k(:,1)+h5*k(:,5)+h6*k(:,6))
call accelerazioni(7,t+h,k)
k(:,10)=y(:,9)
end if
k(:,8) =y(:,9)+h*(bt1*k(:,1)+bt3*k(:,3)+bt4*k(:,4)+bt5*k(:,5))
k(:,9)=y(:,9)+h*(bv1*k(:,1)+bv3*k(:,3)+bv4*k(:,4)+bv5*k(:,5)+bv6*k(:,6))
return
end subroutine rkf
!
subroutine rkfh(k,s,h)
!
! s deve essere compreso tra 0 ed 1 mentre h e' il passo.
!
real(kind=q),intent(in)::s,h
real(kind=q),dimension(1:,1:),intent(out)::k
real(kind=q),parameter::b10=1.0_q,b11=301.0_q/120.0_q,&
b12=-269.0_q/108.0_q,b13=311.0_q/360.0_q
real(kind=q),parameter::b31=7168.0_q/1425.0_q,&
b32=-4096.0_q/513.0_q,b33=14848.0_q/4275.0_q
real(kind=q),parameter::b41=-28561.0_q/8360.0_q, &
b42=199927.0_q/22572.0_q,b43=-371293.0_q/75240.0_q
real(kind=q),parameter::b51=57.0_q/50.0_q,b52=-3.0_q,b53=42.0_q/25.0_q
real(kind=q),parameter::b61=-96.0_q/55.0_q,b62=40.0_q/11.0_q, &
b63=-102.0_q/55.0_q
real(kind=q),parameter::b71=3.0_q/2.0_q,b72=-4.0_q,b73=5.0_q/2.0_q
real(kind=q)::b1,b3,b4,b5,b6,b7
b1=(b10-s*(b11+s*(b12+s*b13)))*s*h
b3=(s*(b31+s*(b32+s*b33)))*s*h
b4=(s*(b41+s*(b42+s*b43)))*s*h
b5=(s*(b51+s*(b52+s*b53)))*s*h
b6=(s*(b61+s*(b62+s*b63)))*s*h
b7=(s*(b71+s*(b72+s*b73)))*s*h
k(:,11)=k(:,10)+b1*k(:,1)+b3*k(:,3)+b4*k(:,4)+b5*k(:,5)+b6*k(:,6)+b7*k(:,7)
return
end subroutine rkfh
!
subroutine accelerazioni_esempio(n,t,k)
!
! Il vettore colonna da utilizzare e' sempre il nono mentre
! il vettore su cui mettere i risultati e' indicato da n.
!
integer,intent(in)::n
real(kind=q),intent(in)::t
real(kind=q),intent(in out),dimension(1:,1:)::k
real(kind=q)::r32
!
if(t>=0.0_q) then
k(1,n)=k(3,8) ! `x
k(2,n)=k(4,8) ! `y
r32=k(1,8)**2+k(2,8)**2
r32=1.0_q/(r32*sqrt(r32))
k(3,n)=-k(1,8)*r32 ! `vx
k(4,n)=-k(2,8)*r32 ! `vy
else
! al tempo 0 valgono queste condizioni iniziali:
! Notare che al tempo zero il parametro del tempo
! serve per dare l'eccentricita' e non il tempo.
! L'inizializzazione dunque vale passando valori negativi
! del tempo. Non e' ammessa una eccentricita' esattamente nulla.
!
r32=min(0.999999_q,abs(t))
k(1,n)=1-r32 ! 1-e
k(2,n)=0.0_q
k(3,n)=0.0_q
k(4,n)=sqrt((1.0_q+r32)/(1.0_q-r32)) ! sqrt((1+e)/(1-e))
write(unit=*,fmt="(1a,4f18.10)") " ",k(:,n)
end if
return
!
! Questo esempio e' il problema del moto di un pianeta attorno al sole
! che occupa l'origine delle coordinate. Il semiasse maggiore e' sempre
! uno, se il moto e' ellittico ossia se l'eccentricita', mai negativa,
! risulta essere inferiore ad 1. Se valesse esattamente 1 il moto
! sarebbe parabolico. Se e indica l'eccentricita', allora a t=0 bisogna
! usare queste condizioni iniziali: k(1,n)=1-e, k(2,n)=0, k(3,n)=0,
! k(4,n)=sqrt((1+e)/(1-e))
!
end subroutine accelerazioni_esempio
!
subroutine faccio_esempio(ec,passi,giri,posizioni)
real(kind=q),intent(in)::ec
integer,intent(in)::passi,giri
real(kind=q),dimension(0:,1:),intent(out)::posizioni
real(kind=q)::t0,s
integer::j,kp
call vedoparametri()
y_esempio(:,:)=0
k_esempio(:,:)=0
t_esempio=0.0_q
h_esempio=2.0_q*atan(1.0_q)/(3.0_q*max(1,passi))
horn_esempio=.true.
call accelerazioni_esempio(9,-abs(ec),y_esempio)
! Ogni j=passi dovrebbe aver fatto un dodicesimo di giro.
do kp=1,giri*12-12
do j=1,passi
call rkf(k_esempio,y_esempio,t_esempio,h_esempio,&
horn_esempio,accelerazioni_esempio)
t_esempio=t_esempio+h_esempio
y_esempio=k_esempio
end do
end do
do kp=1,12
do j=1,passi
call rkf(k_esempio,y_esempio,t_esempio,h_esempio,&
horn_esempio,accelerazioni_esempio)
t_esempio=t_esempio+h_esempio
y_esempio=k_esempio
end do
write(unit=*,fmt="(i5,4f18.10)") kp,y_esempio(:,9)
end do
write(unit=*,fmt=*)"Ora il tempo vale: ",t_esempio
t0=giri*atan(1.0_q)*8.0_q
write(unit=*,fmt=*)"Il tempo voluto sarebbe:",t0
!
do kp=0,630
procedo:do
if(t_esempio>t0) then
exit procedo
end if
! write(unit=*,fmt=*)"Avanzo da",t_esempio," per superare",t0
call rkf(k_esempio,y_esempio,t_esempio,h_esempio,&
horn_esempio,accelerazioni_esempio)
t_esempio=t_esempio+h_esempio
y_esempio=k_esempio
! write(unit=*,fmt="(i5,4f18.10)") kp,y_esempio(:,9)
end do procedo
s=(t0-t_esempio)/h_esempio+1.0_q
! write(unit=*,fmt=*) kp,"Il parametro s=",s
call rkfh(y_esempio,s,h_esempio)
posizioni(kp,1)=y_esempio(1,11)
posizioni(kp,2)=y_esempio(2,11)
! write(unit=*,fmt="(f16.7,4f15.7)") t0,y_esempio(:,11)
! read(unit=*,fmt=*)
t0=t0+0.01_q
end do
return
end subroutine faccio_esempio
!
end module fehlberg
!
!