中心差分法的蛙跳(Leapfrog)格式

差分是用离散形式来近似函数变化率的方法，是数值计算中的一种重要工具。它将函数的连续变化用离散点上的差值表示。中心差分(Central Difference)使用两侧的离散点来近似计算导数。

这就像扁担一样，一前一后，两头都有。当步长为    时

若要求    ，则

式(3)是本文的灵魂所在。

如图1所示，使用蛙跳格式时，速度和加速度、位移不是同一个节拍。根据式(3)，每一个时间步长的中点的速度为

▲图1

式(4)的速度来表示位移

需要指出的是，蛙跳格式以加速度为基本变量，没有直接给出整数时刻的速度。在求解时，先根据运动方程来求解    时刻的加速度

▲图2

以上一篇的小车的自由振动为例，用蛙跳格式模拟。

直观感受动力学数值解法的误差

根据图3，初始条件在    时刻，此时的加速度    可由式(6)求得。    时刻的速度可由向前差分公式(仅用于起步)

在    时刻的位移    可由式(5)求得

接下来便可根据式(6)求    时刻的加速度    ，进一步根据式(4)求    时刻的速度    。求解流程为

▲图3

import numpy as np
import matplotlib.pyplot as plt

# 自由振动蛙跳格式

def main(steps, dt):
    k = 0.5    # 弹簧刚度
    x = 0.5   # 初始位移
    m = 1     #质量
    v1 = 0   # 初始速度
    x1 = x

    xxx = np.zeros(steps)
    time = np.zeros(steps)

    for i in range(steps):
        # 注意蛙跳格式的区别
        a  = -k *x1 /m   # 加速度
        if i == 0 :
            v2 = v1 + a* 0.5* dt   # 起步时用向前差分公式
        else:
            v2 = v1 + a * dt  # 中心差分公式
            
        x2 = x1 + v2*dt # 下一时刻的位移
        
        v1 = v2
        x1 = x2

        # 以下的数据用于绘图
        xxx[i] = x2
        t = dt * i
        time[i] = t

    fig, ax = plt.subplots(1, 1, figsize=(10,6) )
    ax.plot(time, xxx, 'r-')
    ax.set_xlabel('$t/s$',fontsize = 14)
    ax.set_ylabel('$x2/m$',fontsize = 14)

    fig.savefig('f46.png', dpi = 400) 
    plt.show()       



if __name__ == "__main__":
    main(800, 0.1)

▲图4

经过数千次迭代后，小车的振幅仍然接近    。对于该案例而言，蛙跳法提供了良好的模拟效果，但在上万次迭代后误差会变得显著。

蛙跳式中心差分法是显式有限元软件LS-DYNA的最常用的积分方法。