I am trying to build a Runge-Kutta code to integrate the equations of motion for a simple harmonic oscillator. However, when I run the code, I only see first order improvement in the error as I decrease the step size. What could be causing this issue?
Here I start with the second order equation x''+x=0 and introduce a variable v=x'. Then we have a system of two first order equations:
x'=v
v'=-x.
Written in matrix form as (x',v')^T=M(x,v)^T (T stands for transposition), where M=(0,1;-1,0).
Then we compute A,B,C,D in the RK formula x(i+1)=x(i)+h/6*(A+2B+2C+2D) using (ts=timestep)
A=M(x0,v0)
B=M[(x0,v0)+ts/2*A]
C=M[(x0,v0)+ts/2*B]
D=M[(x0,v0)+ts*C]
and apply the RK formula written above. But this gives only linear improvement (ie halving the timestep increases the accuracy by a factor of 2).
program rk4
implicit none
integer :: i,j,k
integer :: nos
double precision, allocatable :: y(:,:)
double precision :: x0,x1,x2,x3,x4,v0,v1,v2,v3,v4
double precision :: ts
allocate(y(2,10000))
x0=1.D0 !initial position
v0=0.D0 !initial velocity
ts=.025D0 !timestep
nos=1000 !number of steps
do i=1,nos
x1=v0
v1=-x0
x2=v0+ts/(2.D0)*v1 !a
v2=-x0-ts/(2.D0)*x1
x3=v0+ts/(2.D0)*v2 !b
v3=-x0-ts/(2.D0)*x2
x4=v0+ts/(2.D0)*v3 !c
v4=-x0-ts/(2.D0)*x3
x0=x0+ts/(6.D0)*(x1+(2.D0)*x2+(2.D0)*x3+x4) !d
v0=v0+ts/(6.D0)*(v1+(2.D0)*v2+(2.D0)*v3+v4)
y(1,i)=x0
y(2,i)=v0
end do
open(1,file='rk4.dat', status='unknown')
do i=1,nos
write(1,*) y(1,i), dble(i)*ts !use to plot a graph of the solution
end do
write(*,*) cos(6.3D0)-y(1,252) !error
end
