! Reshavanje j-ne kretanja za 1D LHO - Runge-Kutta metoda ! Program Program LHO_rk ! zadavanje pochetnih uslova real*8, parameter:: k=1.0d0 real*8, parameter:: x0=1.d0, v0=0.d0, tk=20.d0, dt=0.001, tp=0 !pochetni polozhaj, brzina ! konachno vreme, korak u vremenu, pochetno vreme ! deklaracija tipa svih promenljivih real*8 x(2), y(2), xout(2), energija, yy ! redom : trenutni polozhaj, brzina, ubrzanje, energija, f-ja koja ce se izrachunati real*8 t ! vreme ! Otvaranje fajla (datoteke) u koji ce se smestiti rezultati izrachunavanja open(30,file='lhork4.dat',status='unknown') ! prorachunavanje polozhaja, brzine i energije t=tp x(1)=x0 x(2)=v0 yy=1. ! rachunati sve dok vreme ne bude jednako konachnom vremenu do while(t.le.tk) ! energija lho jedinicne mase energija=1.d0/2.d0*x(2)*x(2)+1.d0/2.d0*k*x(1)*x(1) ! egzaktna trajektorija lho sa parametrima izabranim u knjizi i programu yy=dcos(t) ! ispisivanje podataka u datoteku, 100 je broj naredbe koja pokazuje kako se ispisuju podaci write(30,100) t,x(1),x(2),energija,yy call rk4(x,t,dt) enddo 100 format(10(e13.7,2x)) stop end subroutine derivs(t,x,y) real*8, parameter:: k=1.0d0 real*8 x(2),y(2),t y(1)=x(2) y(2)=-k*x(1) return end ! Runge-Kutta 4tog reda SUBROUTINE rk4(y,x,h) INTEGER, parameter:: n=2,NMAX=2 REAL*8 h,x,dydx(n),y(n),yout(n) !EXTERNAL derivs !Given values for the variables y(1:n) and their derivatives dydx(1:n) known at x, use !the fourth-order Runge-Kutta method to advance the solution over an interval h and return !the incremented variables as yout(1:n), which need not be a distinct array from y. The !user supplies the subroutine derivs(x,y,dydx), which returns derivatives dydx at x. INTEGER i REAL*8 h6,hh,xh,dym(NMAX),dyt(NMAX),yt(NMAX) hh=h*0.5 h6=h/6. call derivs(x,y,dydx) do i=1,n !First step. yt(i)=y(i)+hh*dydx(i) enddo xh=x+hh call derivs(xh,yt,dyt) !Second step. do i=1,n yt(i)=y(i)+hh*dyt(i) enddo call derivs(xh,yt,dym) !Third step. do i=1,n yt(i)=y(i)+h*dym(i) dym(i)=dyt(i)+dym(i) enddo call derivs(x+h,yt,dyt) !Fourth step. do i=1,n !Accumulate increments with proper weights. yout(i)=y(i)+h6*(dydx(i)+dyt(i)+2.*dym(i)) y(i)=yout(i) enddo x=x+h return END