by Reinaldo Baretti Machin
reibaretti2004@yahoo.com
url www.oocities.org/serienumerica
url www.oocities.org/reibaretti2004
 

c Runge Kutta -harmonic oscillator problem -second order equation
c  requires two sets of K functions. The first set develops the
c   position  and we label it K1pos , K2pos etc. The second set
c    develops the velocity. We label it M1vel , M2vel etc.
cThe chosen values for k and m give omega=1.rad/s . The period is 2 pi.
c Initial conditions are x1=1. and velocity (x2=0.). The analytic
c  solution is x=cos(omega*t)  , v=-omega*sin(omega*t)
c     x1 stands for position , x2 represents the velocity

      real k,m ,k1pos,k2pos,k3pos,k4pos,m1vel,m2vel,m3vel,m4vel
      f2(t,x1,x2)= -(k/m)*x1
      f1(t,x1,x2)= x2
      pi=2.*asin(1.)
      x1=1.
      x2=0.
      pos0=x1
      vel0=x2
      k=1.
      m=1.
      omega=sqrt(k/m)
      period=2.*Pi/omega
      kount=10
      h=period/(100.)
      do 10 i =1,100
      t=h*float(i)
      k1pos=h*f1(t,x1,x2)
      m1vel=h*f2(t,x1,x2)
c*******************
      k2pos=h*f1(t+h/2.,x1+k1pos/2.,x2+m1vel/2.)
      m2vel=h*f2(t+h/2.,x1+k1pos/2.,x2+m1vel/2.)
c*****************************
      k3pos=h*f1(t+h/2., x1+k2pos/2., x2+m2vel/2.)
      m3vel=h*f2(t+h/2., x1+k2pos/2., x2+m2vel/2.)
c***************************
      k4pos=h*f1(t+h,x1+k3pos,x2+m3vel)
      m4vel=h*f2(t+h,x1+k3pos,x2+m3vel)
c*************************
      pos1=pos0+(1./6.)*(k1pos+2.*k2pos+2.*k3pos+k4pos)
      vel1=vel0+ (1./6.)*(m1vel+2.*m2vel+2.*m3vel+m4vel)
      If(i.eq.kount)then
      
      print*,'t, pos,vel=',  t, pos1, vel1
      print*,'a cos(wt),-w*asin(wt)=',cos(omega*t),-omega*sin(omega*t)
      kount=kount+10
      print*,'   '
      endif
      x1=pos1
      x2=vel1
      pos0=pos1
      vel0=vel1
10    continue
      stop
      end
      
      RUN
      
       t, pos,vel=  0.628318548  0.809017062 -0.587785244
 a cos(wt),-w*asin(wt)=  0.809017003 -0.587785244

 t, pos,vel=  1.2566371  0.309017062 -0.95105654
 a cos(wt),-w*asin(wt)=  0.309016973 -0.95105654

 t, pos,vel=  1.88495564 -0.309016913 -0.9510566
 a cos(wt),-w*asin(wt)= -0.309017032 -0.95105648

 t, pos,vel=  2.51327419 -0.809017003 -0.587785363
 a cos(wt),-w*asin(wt)= -0.809017062 -0.587785184

 t, pos,vel=  3.14159274 -1.00000012 -9.00191921E-008
 a cos(wt),-w*asin(wt)= -1.  8.74227766E-008

 t, pos,vel=  3.76991129 -0.809017241  0.587785244
 a cos(wt),-w*asin(wt)= -0.809016943  0.587785363

 t, pos,vel=  4.39823008 -0.309017211  0.951056659
 a cos(wt),-w*asin(wt)= -0.309016645  0.9510566

 t, pos,vel=  5.02654839  0.309016883  0.951056838
 a cos(wt),-w*asin(wt)=  0.309017122  0.95105648

 t, pos,vel=  5.65486717  0.809017062  0.587785602
 a cos(wt),-w*asin(wt)=  0.809017241  0.587784946

 t, pos,vel=  6.28318548  1.00000036  2.83734295E-007
 a cos(wt),-w*asin(wt)=  1. -1.74845553E-007