Haskell中的Verlet集成

Question

我是Haskell的新手,作为一个练习,我一直在努力实现Joel Franklin的" 计算方法在物理学"一书中的一些代码(用Mathematica编写).我编写了以下代码,将lambda表达式(加速度)作为第一个参数.通常,加速度的形式为x''= f(x',x,t),但并不总是所有三个变量.

-- Implementation 2.1: The Verlet Method
verlet _  _  _  _  0 = [] 
verlet ac x0 v0 dt n = take n $ zip [0,dt..] $ verlet' ac x0 v0 0 dt
  where verlet' ac x0 v0 t0 dt = 
          let
            xt = x0 + dt*v0 + 0.5*(dt^2)*(ac x0 v0 t0)
            vt = v0 + dt*(ac x0 v0 t0)
          in
            xt:(verlet' ac xt vt (t0+dt) dt)

在ghci中,我将使用以下命令运行此代码(加速函数a = - (2pi)² x来自书中的练习):

verlet (\x v t -> -(2*pi)^2*x) 1 0 0.01 1000

我的问题是这不是真正的Verlet方法 - 这里x _{n + 1} = x _n + dt*v _n + 1/2*a(x _n,v _n,n),而维基百科中描述的Verlet方法去x _{n + 1} = 2*x _n - x _n-1 + a(x _n,v _n,n).如何重新编写此函数以更忠实地表示Verlet集成方法？

从切边来看,有没有办法更优雅,更简洁地写出这个？是否有线性代数库可以使这更容易？我感谢你的建议.

Answer 1

忠实的Verlet序列具有x _n,取决于x-x _n-1和x _n-2的前两个值.这种序列的典型例子是Fibonacci序列,它具有这种单行Haskell定义:

fibs :: [Int]
fibs = 0 : 1 : zipWith (+) fibs     (tail fibs)
                        -- f_(n-1)  f_n

这将Fibonacci序列定义为无限(懒惰)列表.自我引用tail fibs为您提供前一个术语,并且引用fibs为您提供之前的术语.然后将这些术语组合(+)以产生序列中的下一个术语.

您可以采用与您的问题相同的方法,如下所示:

type State = (Double, Double, Double)  -- (x, v, t) -- whatever you need here

step :: State -> State -> State
step s0 s1 = -- determine the next state based on the previous two states

verlet :: State -> State -> [State]
verlet s0 s1 = ss
  where ss = s0 : s1 : zipWith step ss (tail ss)

数据结构State包含您需要的任何状态变量 - x,v,t,n,......该函数step类似于(+)Fibonacci情况,并计算给定前两个的下一个状态.verlet给定初始两个状态时,该函数确定整个状态序列.