Do wyznaczenia polozen obietku w dowolnej chwili czasu jest metoda Rungego-Kutty. To metoda dokladniejsza niz poprzednia przedstwawiona, metoda Eulera. Opiera sie na rozwinieciu wyrazow w szereg do czwartego rzedu. Ponizszy kod programu ORBIT1 przedstawia te metode:
program orbit1
implicit none
real epsilon,ykonc(4)
real y0(4),y1(4),y2(4),y3(4),f(4),f1(4),f2(4),f3(4)
real k1(4),k2(4),k3(4),k4(4)
integer i, j,kroki,mn
open(22,file='orbita1.dat', status='unknown')
write(
,)'Podaj ilosc krokow:'
read(
,) kroki
epsilon=0.01 podaje krok czasowy
y0(1)=0.5 podaje wartosc x0
y0(2)=0.0 podaje wartosc y0
y0(3)=0.0 podaje wartosc Vx
y0(4)=1.63 podaje wartosc Vy
write (22,*)y0(1), y0(2)
wywyoluje procedury:
do j=1,kroki
call jeden(y0,f) oblicza funkcje F1
call dwa(f,k1,epsilon) oblicza wsp/o/lczynnik K1
call trzy(k1,0.5,y0,y1) oblicza y1
call jeden(y1,f1) oblicza funkje F2
call dwa(f1,k2,epsilon) oblicza wsp/o/lczynnik K2
call trzy(k2,0.5,y0,y2) oblicza y2
call jeden(y2,f2) oblicza funkcje F3
call dwa(f2,k3,epsilon) oblicza wsp/o/lczynnik K3
call trzy(k3,1.0,y0,y3) oblicza y3
call jeden(y3,f3) oblicza funkcje F4
call dwa(f3,k4,epsilon) oblicza wspolczynnik K4
do i=1,4
ykonc(i)=y0(i)+(1./6.)*k1(i)+(1./3.)*k2(i)+ oblicza ostateczn/a
& (1./3.)*k3(i)+(1./6.)*k4(i) warto/s/c wektora Y
enddo
write(22,*)ykonc(1), ykonc(2)
do i=1,4
y0(i)=ykonc(i)
enddo
enddo
end
c +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
subroutine jeden(y,f)
c liczy funkcje F(1...4)
real y(4), f(4),r
f(1)=y(3)
f(2)=y(4)
r=sqrt(y(1)**2.0+y(2)**2.0)
f(3)=-y(1)/r**3.0
f(4)=-y(2)/r**3.0
return
end
subroutine dwa(f,k,epsilon)
c oblicza wsolczynniki K
real k(4), f(4), epsilon
k(1)=epsilon*f(1)
k(2)=epsilon*f(2)
k(3)=epsilon*f(3)
k(4)=epsilon*f(4)
return
end
subroutine trzy(k,wspol,y0,y1)
c oblicza y(1...4)
real k(4),wspol,y0(4),y1(4)
y1(1)=y0(1)+wspol*k(1)
y1(2)=y0(2)+wspol*k(2)
y1(3)=y0(3)+wspol*k(3)
y1(4)=y0(4)+wspol*k(4)
return
end