2
$\begingroup$

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
$\endgroup$
7
$\begingroup$

It appears you are only taking half a time step on the fourth stage (where you should take a full one). Since you aren't taking the proper time step, this is turning into a complicated but relatively useless explicit scheme.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.