5

我是 Haskell 的新手，作为练习，我一直在尝试实现 Joel Franklin 的著作Computational Methods in Physics中的一些代码（用 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) 2 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 集成方法？

切线，有没有办法更优雅和简洁地写这个？是否有线性代数库可以使这更容易？我很感激你的建议。

4

3 回答 3

5

忠实的 Verlet 序列具有 x n取决于 x 的前两个值 - x n-1和 x n-2。这种序列的一个典型例子是斐波那契序列，它具有这个单线 Haskell 定义：

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

这将斐波那契数列定义为无限（惰性）列表。自引用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类似于(+)斐波那契情况，并根据前两个状态计算下一个状态。该verlet函数确定给定初始两个状态的整个状态序列。

于 2014-03-09T06:50:12.770 回答
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) 是某个势函数的梯度，而势函数是静态对象，它们只依赖于位置，而不是时间和速度。许多无摩擦机械系统都属于这一类。

速度维莱特

现在使用一阶中心差商将时间 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

这是编写velocity-Verlet 算法的一种变体。

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

使用上一步从 n-1 到 n 的方程，可以在速度计算中将 x n-1替换为 v n-1和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进行碰撞检测

一个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 0 = a( x 0 )。

跳蛙Verlet

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

一个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

一个n = a(x n )

v n+0.5 = v n-0.5 + a n *dt

初始化只需要计算

v 0+0.5 = v 0 + 0.5*a( x 0 )*dt

（结束更新）

由于时间索引不匹配，使用 x n和 v n-0.5或 v n+0.5计算的任何统计数据都会因与 dt 成比例的误差而偏离。碰撞可以单独使用位置矢量来检测，但在偏转中，速度也需要进行明智的更新。

于 2014-03-09T17:50:08.993 回答
0

这是我在实施 user5402 的建议后的解决方案：

-- 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 中做到这一点有任何建议，欢迎他们）。此外，如果您看到任何明显的重构，请随时指出它们。谢谢！

于 2014-03-10T07:04:43.130 回答