Fortran 中 Verlet 算法的错误顺序错误

Question

我正在尝试使用 Fortran 中的 Verlet 方法（Original Verlet）来模拟谐波振荡器。我的研究表明错误的阶数应该是 2，但我的计算显示是 1 的阶数。

我在我的源代码中找不到我的错误。我该怎么办？

编辑：我使用的算法如下：

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

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

其中x(t)代表位置，v(t)代表速度，F(t)代表力。我认出这是描述的原始 Verlet here

根据this site，错误的顺序至少应该是O(Δt²) 但是我在gnuplot 中绘制的程序顺序错误（下图）没有O(Δt²) 的顺序.

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

以上程序计算了时间步长h的计算误差。

Answer 1

根据 Wiki page for Verlet integrators，似乎我们需要使用更准确的方法来设置 xold 的初始值（即包括力项）以得到 2 阶的全局误差。确实，如果我们将 xold 修改为

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

全局误差变为 2 阶（请参见下面的对数-对数图）。

Fortran 中 Verlet 算法的错误顺序错误

Incorrect Order of Error of Verlet Algorithm in Fortran

fortran

numerical-methods