我是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)
verlet (x v t -> -(2*pi)^2*x) 1 0 0.01 1000
我的问题是，这不是真正的Verlet方法-这里xn+1 = xn + dt*vn+1/2*a(xn，vn，n)，而维基百科中描述的Verlet方法是xn+1 = 2*xn - xn，vn，n)。我该如何重写这个函数，以更忠实地表示Verlet集成方法?
顺便问一句，有没有一种更优雅、更简洁的方式来写这个?有没有线性代数库能让这更简单?谢谢你的建议。忠实Verlet序列有xn取决于x的前两个值——xn-1和xn-2。这种序列的一个典型例子是Fibonacci序列，它有这样一行的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类似于Fibonacci中的(+)，并根据前两个状态计算下一个状态。verlet函数确定给定初始两个状态的整个状态序列。
实际上，如果你继续读下去，你会发现这两种变体在维基百科页面上都有。
基本Verlet
二阶ODE x''(t)=a(x(t))的基本二阶中心差商离散为
x <子> n + 1 n - 2 * x <子> + x <子> n - 1 = <子> n * dt ^ 2
注意，迭代中没有速度，加速度函数a(x)中也没有速度。这是因为Verlet积分只有在动力系统是保守的情况下才优于其他积分方法，即-m*a(x)是某个势函数的梯度，而势函数是静态对象，它们只与位置有关，与时间和速度无关。许多无摩擦机械系统都属于这一类。
速度Verlet
现在，利用一阶中心差商，设tn时刻的速度为
v <子> n * (2 * dt) = x <子> n + 1 - x <子> n - 1
与第一个方程加减，得到
2 * x <子> n + 2 * x <子> n - 1 = 2 * v <子> n * dt + n <子> * dt ^ 2
2 * x <子> n + 1 n - 2 * x <子> = 2 * v <子> n * dt + n <子> * dt ^ 2
或
v <子> n = (x <子> n - x <子> n - 1 )/dt + 0.5 * n <子> * dt
x <子> n + 1子> n n + v <子> * dt + 0.5 * n <子> * dt ^ 2
这是velocity-Verlet算法的一种变体。
(更新为使所有状态变量在迭代步骤之前和之后对应于相同的时间)
利用前一步从n-1到n的方程，可以用vn-1和an-1来代替xn-1。
v <子> n = v <子> n - 1 + 0.5 *(<子> n - 1 + n <子> ) * dt
为了避免任何向量x,v,a的两个实例，可以安排更新过程，使所有内容都到位。假设在迭代步骤开始时，存储的数据对应于(tn-1，xn-1，vn-1，an-1)。然后计算下一个状态为
v <子> n - 0.5 v = <子> n - 1 + 0.5 * <子> n - 1子> p> x <子> n = x <子> n - 1 + v <子> n - 0.5 * dt
使用xn和vn-0.5
进行碰撞检测an = a(xn)
v <子> n = v <子> n - 0.5 + 0.5 * n <子> * dt
使用xn和vn
进行统计
或作为代码
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方案并显著改变结果，操作的旋转是可能的，但必须小心迭代步骤后对状态的解释。
唯一需要初始化的是计算a0 = a(x0)。超越Verlet
从速度Verlet公式可以看出，对于位置的更新，不需要速度vn，而只需要半点速度vn+0.5。
an = a(xn)
0.5 v <子> n + v = <子> n - 0.5 + n <子> * dt
x <子> n + 1子> n + v <子> n + 0.5 * dt
或代码
a = eval_a(x);
v += a*dt;
x += v*dt;
同样，这些操作的顺序是非常重要的，对于保守系统，改变顺序会导致奇怪的结果。
(Update)然而，可以将执行顺序旋转到
x += v*dt;
a = eval_a(x);
v += a*dt;
这对应于三元组(tn，xn，vn+0.5)的迭代，如
x <子> n = x <子> n - 1 + v <子> n - 0.5 * dt
an = a(xn)
0.5 v <子> n + v = <子> n - 0.5 + n <子> * dt
初始化只需计算
v <子> 0 + 0.5 v = <子> 0 + 0.5 * (x <子> 0 )* dt
(结束更新)
使用xn和vn-0.5或vn+0.5计算的任何统计数据都将由于与dt成正比的误差而关闭，因为时间指标不匹配。碰撞可以用位置矢量单独检测，但在偏转中，速度也需要合理地更新。
以下是我在执行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中做到这一点，欢迎他们)。此外，如果您看到任何明显的重构，请随时指出它们。谢谢!