Haskell中的Verlet集成

忠实的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 函数确定给定初始两个状态的整个状态序列。

实际上,如果你继续阅读,你会发现这两种变体都在维基百科页面上显示。

基本Verlet

二阶ODE x“”(t)=a(x(t))的基本二阶中心差商离散化为

x(x) _n+1个 -2*x个 _n +x(x) _n-1个 =一个 _n *时间^2

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

Verlet速度

现在,使用一阶中心差商,设置时间t的速度 _n 到

v _n *(2*dt)=x _n+1个 -x(x) _n-1个

并将其与第一个方程相加和相减,以获得

-2*x个 _n +2*x个 _n-1个 =-2*伏 _n *时间+a _n *时间^2

2*x个 _n+1个 -2*x个 _n =2*v _n *时间+a _n *时间^2

或

v _n =(x _n -x(x) _n-1个 )/dt+0.5*a _n *时间

x(x) _n+1个 =x _n +五 _n *dt+0.5*a _n *时间^2

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

( 更新以使所有状态变量对应于迭代步骤前后的相同时间 )

使用前一步从n-1到n的方程式,可以替换x _n-1个 通过v _n-1个 和一个 _n-1个 在速度计算中。然后

v _n =伏 _n-1个 +0.5*(a _n-1个 +一个 _n )*时间

为了避免任何向量x、v、a都有两个实例,可以安排更新过程,以便一切就绪。假设在迭代步骤的入口,存储的数据对应于(t _n-1个 ,x _n-1个 v _n-1个 一 _n-1个 ). 然后下一个状态被计算为

v _正0.5 =伏 _n-1个 +0.5*年 _n-1个 *时间

x(x) _n =x _n-1个 +五 _正0.5 *时间

使用x进行碰撞检测 _n 和v _正0.5

一 _n =a(x _n )

v _n =伏 _正0.5 +0.5*年 _n *时间

使用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(x ₀ ).

Leapfrog Verlet公司

从速度Verlet的公式可以看出,对于位置的更新,不需要速度v _n 但只有半点速度v _n+0.5 然后

一 _n =a(x _n )

v _n+0.5 =伏 _正0.5 +一个 _n *时间

x(x) _n+1个 =x _n +五 _n+0.5 *时间

或在代码中

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(x) _n =x _n-1个 +五 _正0.5 *时间

一 _n =a(x _n )

v _n+0.5 =伏 _正0.5 +一个 _n *时间

初始化只需要计算

v _0+0.5 =伏 ₀ +0.5*a(x ₀ )*时间

( 结束更新 )

使用x计算的任何统计信息 _n 和v _正0.5 或v _n+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中绘制x有任何建议,欢迎他们)。此外,如果您看到任何明显的重构,请随时指出它们。谢谢