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