1

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.

Lutz Lehmann
  • 25,219
  • 2
  • 22
  • 51
Kaira
  • 133
  • 4
  • 1
    The Verlet method only works with all the claimed properties if the force is a function of `x`, and furthermore is the negative gradient of a potential energy function. Meaning it should be `F(x(t))`. – Lutz Lehmann Jul 28 '19 at 08:45
  • I am curious to know the need for `real128` (real*16) as opposed to double-precision `real64`. That alone would likely more than double the runtime of your code. Also, `real*16` is not Fortran-standard conforming. A formal way would be to `use iso_fortran_env, only: real128` and declare real variables as `real128`. – Scientist Jul 28 '19 at 17:58

1 Answers1

2

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

roygvib
  • 7,218
  • 2
  • 19
  • 36
  • @IanBush Yes, I think you are right.. Though I tried to emphasize the key modified part (while not changing other parts), I'll add a link for OP below. – roygvib Jul 29 '19 at 08:27
  • @MaoNishino As mentioned above, real*16 is not portable, so it is recommended to use the KIND mechanism; plsease see, for example, this Q/A as a starter: https://stackoverflow.com/questions/3170239/fortran-integer4-vs-integer4-vs-integerkind-4 Also, please note that "1.23" etc means default real type (typically 32 bits). We can specify the kind of such a literal as 1.23_dp etc (where dp is a KIND parameter). – roygvib Jul 29 '19 at 08:29