Search code examples
fortrannumerical-methods

Incorrect Order of Error of Verlet Algorithm in Fortran

I am trying to simulate harmonic oscillator by using Verlet Method(Original Verlet) in Fortran. My research tells that the order of error should be 2 but my calculation showed the order of 1.

I couldn't find my mistake in my source code. What should I do?

Edit: The algorithm I am using is below:

x(t+Δt) = 2x(t) - x(t-Δt) + Δt² F(t)/m

v(t) = {x(t+Δt) -x(t-Δt)}/2Δt

Where x(t) represents the position, v(t) represents velocity and F(t) represents Force. I recognize this is the Original Verlet described here

According to this site, the order of error should be at least O(Δt²) but the error of the order of my program plotted in gnuplot (below) does not have a order of O(Δt²).

plot

program newton_verlet

implicit none

real*16, parameter :: DT   = 3.0
real*16, parameter :: T0   = 0.0
real*16, parameter :: TEND = 2.0
integer, parameter :: NT   = int(TEND/DT + 0.5)

real*16, parameter :: M    = 1.0
real*16, parameter :: X0   = 1.0
real*16, parameter :: V0   = 0.0

real*16 x,v,t,xold,xnew,vnew,ek,ep,et,f,h
integer it,n

do n=0,20
h = DT/2**n

x = X0
v = V0

ek = 0.5*M*v*v
ep = x*x/2
et = ek + ep

xold = x - v*h

do it = 1,2**n
!   f = -x
   f = -x

   xnew = 2.0*x - xold + f*h*h/M
   v = (xnew-xold)/(2.0*h)

   ek = 0.5*M*v*v
   ep = 0.5*xnew*xnew
   et = ek + ep
   xold = x
   x = xnew

end do
write(*,*) h,abs(x-cos(DT))+abs(v+sin(DT))
end do

end program

Above program calculates the error of calculation for the time step h.

Solution

  • According to the Wiki page for Verlet integrators, it seems that we need to use a more accurate way of setting the initial value of xold (i.e., include terms up to the force) to get the global error of order 2. Indeed, if we modify xold as

    program newton_verlet
implicit none
real*16, parameter :: DT   = 3.0
real*16, parameter :: M    = 1.0
real*16, parameter :: X0   = 1.0
real*16, parameter :: V0   = 0.0

real*16 x,v,xold,xnew,f,h
integer it,n

do n = 0, 20

    h = DT / 2**n

    x = X0
    v = V0
    f = -x

    ! xold = x - v * h                     !! original
    xold = x - v * h + f * h**2 / (2 * M)  !! modified

    do it = 1, 2**n

        f = -x

        xnew = 2 * x - xold + f * h * h / M

        xold = x
        x = xnew
    end do

    write(*,*) log10( h ), log10( abs(x - cos(DT)) )
end do

end program

    the global error becomes of order 2 (please see the log-log plot below).

    compare.png