Verlet Integration在Haskell中的应用

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积分方法？

另外，有没有一种更优雅简洁的编写方法？是否有线性代数库可以使这更容易？谢谢您的建议。

- castle-bravo

Answer 1

忠实的Verlet序列由x_n-1和x_n-2之前两个值决定x_n。其中一个经典例子是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根据初始两个状态确定整个状态序列。

- ErikR

---

这是一个不错的解决方案，但当我尝试将其加载到ghci中进行测试时，我收到了以下错误信息：在数据/新类型声明中的构造函数中出现解析错误：（Double，Double，Double）我会将此帖移到另一个线程。 - castle-bravo

@castle-bravo：答案中有一个错别字（我已经修复了）。user5402在应该使用“type”时使用了“data”；后者用于定义类型同义词（旧类型的新名称），而“data”用于定义全新的类型，因此需要构造函数名称。（这将是另一个选项：“data State = State Double Double Double”，或“data State = State {position，velocity，time :: Double}”。）如果您想要了解更多信息，有一个关于定义类型的Learn You A Haskell章节，如果有任何不清楚的地方。 - Antal Spector-Zabusky

@AntalS-Z和user5402 - 感谢你们的帮助！我已经在下面的回答中放置了你们建议的实现。 - castle-bravo

Answer 2

实际上，如果你继续阅读，你会发现维基百科页面中两种变体都有所涉及。

---

基本Verlet算法

对于二阶ODE x''(t)=a(x(t))的基本二阶中心差商离散化为：

x_n+1 - 2*x_n + x_n-1 = a_n*dt^2

请注意，在迭代中没有速度，也没有加速度函数a(x)。这是因为当动力系统是保守时，Verlet积分仅优于其他积分方法，即-m*a(x)是某个势函数的梯度，并且势函数是静态对象，它们仅取决于位置，而不取决于时间和速度。许多无摩擦机械系统属于此类别。

---

速度Verlet算法

现在使用一阶中心差商，将时间t_n处的速度设置为

v_n*(2*dt) = x_n+1 - x_n-1

并将其添加到第一个方程中，并从中减去，以获得

-2*x_n + 2*x_n-1 = -2*v_n*dt + a_n*dt^2

2*x_n+1 - 2*x_n = 2*v_n*dt + a_n*dt^2

或者

v_n = (x_n - x_n-1)/dt + 0.5*a_n*dt

x_n+1 = x_n + v_n*dt + 0.5*a_n*dt^2

这是一种编写速度Verlet算法的变体。

---

（更新为在迭代步骤之前和之后所有状态变量对应于相同的时间）

使用上一步的方程从n-1到n，可以将x_n-1替换为v_n-1和a_n-1以计算速度。然后：

v_n = v_n-1 + 0.5*(a_n-1 + a_n)*dt

为避免向量x、v、a的任意两个实例，可以安排更新过程使一切就位。假设在迭代步骤的输入时，存储的数据对应于(t_n-1,x_n-1,v_n-1,a_n-1)。然后计算下一个状态为：

v_n-0.5 = v_n-1 + 0.5*a_n-1*dt x_n = x_n-1 + v_n-0.5*dt 使用x_n和v_n-0.5进行碰撞检测 a_n = a(x_n) v_n = v_n-0.5 + 0.5*a_n*dt 使用x_n和v_n进行统计

或者作为代码：

v += a*0.5*dt;
x += v*dt;
do_collisions(x,v);
a = eval_a(x);
v += a*0.5*dt;
do_statistics(x,v); 

改变这些操作的顺序将会破坏Verlet方案并显着改变结果，可以旋转操作，但必须小心迭代步骤后状态的解释。

唯一需要初始化的是计算a₀=a(x₀)。

跳跃Verlet

从速度Verlet公式中可以看出，对于位置的更新，只需半点速度v_n+0.5，而不需要速度v_n。然后：

a_n=a(x_n) v_n+0.5=v_n-0.5+a_n*dt x_n+1=x_n+v_n+0.5*dt

或者在代码中：

a = eval_a(x);
v += a*dt;
x += v*dt;

这些操作的顺序非常重要，对于保守系统来说，更改会导致奇怪的结果。

---

(更新) 然而，可以旋转执行顺序以

x += v*dt;
a = eval_a(x);
v += a*dt;

这对应于三元组(t_n,x_n,v_n+0.5)的迭代，如下所示

x_n = x_n-1 + v_n-0.5*dt

a_n = a(x_n)

v_n+0.5 = v_n-0.5 + a_n*dt

初始化只需要计算

v_0+0.5 = v₀ + 0.5*a( x₀ )*dt

(结束更新)

使用x_n和v_n-0.5或v_n+0.5计算的任何统计量都会存在误差，其误差与时间步长dt成正比，因为时间索引不匹配。位置向量可以检测到碰撞，但在偏转时速度也需要得到合理的更新。

- Lutz Lehmann

Answer 3

在参考了用户5402的建议后，这是我的解决方案：

-- 1-Dimensional Verlet Method
type State = (,,) Double Double Double -- x, x', t

first :: State -> Double
first (x, _, _) = x

second :: State -> Double
second (_, x, _) = x

third :: State -> Double
third (_, _, x) = x

verlet1 :: (State -> Double) -> State -> Double -> Int -> [State]
verlet1 f s0 dt n = take n ss
  where
    ss = s0 : s1 : zipWith step ss (tail ss)
      where 
        s1 = (first s0 + dt*(second s0) + 0.5*dt^2*(f s0), 
              second s0 + dt*(f s0), third s0 + dt)
        step :: State -> State -> State
        step s0 s1 = (2*(first s1) - first s0 + dt^2*(f s1), 
                      second s1 + dt*(f s1), third s1 + dt)

我使用以下命令在ghci中运行它：

verlet1 (\x -> -(2*pi)^2*(first x)) (1, 0, 0) 0.01 100

这似乎产生了我期望的结果-显然是正弦运动！ 我还没有绘制x（如果有人知道如何在Haskell中做到这一点，欢迎提出建议）。 另外，如果您看到任何明显的重构，请随时指出。 谢谢！

- castle-bravo

---

这似乎是速度和跳跃Verlet的奇怪混合。从第一步可以看出，每个状态s = (t, x, v)对应于s[n] = (t[n], x[n], v[n+0.5])。然后，s[0]的初始化应该为x[0]=x0, t[0]=t0, v[0+0.5] = v0 + 0.5*dt*(f (x0, v0, t0))。利用这个，步骤可以简化为跳跃方式：step s = ( first s + dt*(second s), second s + dt*(f s), third s + dt)。 - Lutz Lehmann

或者我读错了，状态应该是s[n]=(t[n], x[n], v[n-0.5])。那么如果想要避免对f的冗余调用，则需要先更新v[n+0.5] = v[n]+dt*(f s[n])，然后按照完全相同的顺序更新x[n+1] = x[n]+dt*v[n+0.5]。 忽略关于建立状态的重复部分，可以仅使用上一个状态来编写步骤，如：step s = ( first s + dt*(second s + dt*(f s)), second s + dt*(f s), third s + dt)。 - Lutz Lehmann