1 常微分方程及其数值解法1.1 常微分方程概述在数学上，物质的运动和变化规律是通过函数关系来表示的，在一些复杂的现象中，我们要求的未知量就变成了满足特定条件的一个或一些未知函数。

有的时候，我们需要利用导数或者微分的关系，即这些未知函数的导数与自变量满足某种关系，这种方程我们称之为微分方程。

未知函数是一元函数的微分方程称之为常微分方程，未知函数是多元函数的微分方程我们称之为偏微分方程，我们这里只考虑常微分方程。

常微分方程的解，就是找出一个代入方程使之成为恒等式的函数。

若微分方程的解中含有任意常数的个数与方程阶数相同，且任意常数之间不能合并，则称此解为该方程的通解。

当通解中的各任意常数都取特定值时所得到的解，称为方程的特解。

在实际问题中，这些函数往往还需要满足一些特定条件，这称之为定解条件。

但在实际问题中，很多常微分方程的解析表达式过于复杂，甚至得不到通解的解析表达式。

而且，常微分方程的特解是否存在，存在几个特解，这涉及到微分方程解的存在性和唯一性定理。

因此，在实际应用中，我们通常利用数值的方法来求得方程的数值解，在误差允许的范围内，我们用数值解来替代解析解。

所以，研究常微分的数值解法是很有必要的。

2.2 常微分方程的数值解法常微分方程的数值解法是有常微分方程的定解条件提出的，首先我们考虑如下一阶常微分方程的初值问题。

()()00(,)dx t f x t dtx t x⎧=⎪⎨⎪=⎩（2.1） 2.2.1 欧拉法欧拉法（又称差分法）是常微分方程初值问题数值解法中最简单最古老的方法，它的基本思路是将（2.1）式中导数项用差分来逼近，从而将一个微分方程转化为一个代数方程，以便迭代求解。

根据用于逼近的差分方式来分，可以分为向前差分、向后差分、中心差分。

()()()()()()()()()111112l l l l l l l l l dx t x t x t dt tdx t x t x t dt tdx t x t x t dt t++++--=∆-=∆-=∆ （2.2） 上式中，分别为向前差分法、向后差分法、中心差分法。

以向前差分法为例，我们在这里推导向前欧拉法的迭代格式。

由（2.2）式中第一个向前差分法公式我们可以得出()()()1l l l dx t x t x t tdt+=+∆ （2.3） 再由原问题（2.1）可以得到向前欧拉法的迭代格式为：()()()1,l l l x t x t tf x t +=+∆ （2.4）因此我们可以利用初始条件，依次得到后一个时刻方程的解，直至求得的解收敛到所要求的时刻为止。

当原问题为更为复杂的常微分方程组或者高阶常微分时，只需要将x 看做向量，原问题就可以简化或者降阶为一个一阶常微分方程组或者一个一阶常微分方程。

这样就可以利用欧拉法的迭代格式进行迭代计算。

2.2.2 改进欧拉法用数值积分的方法对原问题（2.1）式进行离散化处理，方程两边做积分有：()()()11(),l lt l l t x t x t f x t t dt ++-=⎰（2.5）对右端积分使用梯形积分公式可得()()()()()111(),,,2l lt l l l l t tf x t t dt f x t t f x t t +++∆⎡⎤≈+⎣⎦⎰（2.6） 将（2.6）式带入到（2.5）式中可以得到()()()()()()111,,2l l l l l l tx t x t f x t t f x t t +++∆⎡⎤=++⎣⎦ （2.7） 在实际计算过程中，改进的欧拉法是用欧拉法先求得一个预测值()1l xt +，再利用这个预测值来来计算()1l x t +，因此改进欧拉法的迭代格式为：()()()()()()()()()()1111,,,2l l l l l l l l l l x t x t tf x t t tx t x t f x t t f x t t ++++⎧=+∆⎪⎨∆⎡⎤=++⎪⎣⎦⎩ （2.8） 从（2.7）式不难看出改进欧拉法是隐式迭代算法，而欧拉法是显示迭代算法。

2.2.3 Newmark-beta 算法Newmark-beta 算法是用于解决微分方程的数值积分算法之一，它通常用于有限单元法中的动力学分析。

根据运动学方程：212x xt xt =+（2.7） 利用广义中值定理，Newmark-beta 算法对于一次微分项修改为()()()1l l x t x t tx t γ+=+∆ （2.8）和欧拉法不同的是，Newmark-beta 算法对于二次微分项不是直接利用上一个迭代步中的二次微分结果，而是对当前迭代步和上一迭代步两个迭代步的二次微分结果进行了加权取值，利用加权修正后的二次微分结果进行迭代计算。

()()()11() 01l l x t x t x t γγγγ+=-+≤≤ （2.9）将（2.9）式带入到（2.8）式中得到，一次微分项的迭代格式为：()()()()()111l l l l x t x t tx t tx t γγ++=+-∆+∆ （2.10）同理，对原函数也进行相同的处理，利用广义中值定理对位移进行修正得到：()()()()2112l l l x t x t tx t t x t β+=+∆+∆ （2.11）依旧利用加权来修正上式中的二次微分项为()()()()1122 01l l x t x t x t ββββ+=-+≤≤ （2.12）现假设Newmark-beta 算法中令0,1/2βγ==得到，Newmark-beta 算法的迭代格式为()()()()()()()()()()()()()111211222l l l l l l l l l l x t x t x t tx t x t x t x t t x t x t tx t x t γ++++=+∆=++∆=+∆+ （2.13）由上式可以利用初始条件依次得到每个迭代步的数值解。

2.2.4 Leap-Frog （蛙跳）算法Leap-Frog （蛙跳）算法是对一次微分项即速度进行修正，它的基本原理是，首先利用当前时刻的加速度，计算半个时间步长后的速度：()()()1/21/2l l l x t x t x t t +-=+∆ （2.14）之后，计算下一个步长时刻的位置：()()()11/2l l l x t x t x t t ++=+∆ （2.15）利用两个半个时间步长的速度来计算当前时刻的速度()()()1/21/22l l l x t x t x t -++=（2.16）在迭代启动的时候需要对初始条件进行修改，需要用初始条件来反算()1/2xt -：()()()1/200/2x t x t x t t -=-∆ （2.17）2.2.5 Velocity Verlet 算法Velocity Verlet 算法由于其速度计算精度更高被广泛用在分子动力学计算中，它的基本原理为，对位移进Taylor 展开，其展开式为：()()()()()2311126l l l l i l x t x t x t t x t t b t t +=+∆+∆+∆ （2.18） 略去高阶项就得到了Velocity Verlet 算法的迭代格式()()()()2112l l l l x t x t x t t x t t +=+∆+∆ （2.19a ） ()()()()1112l l l l x t x t x t x t t ++=++∆⎡⎤⎣⎦ （2.19b ） 还等价于()()()()()()()()()1/211/211/211212l l l l l l l l l x t x t x t tx t x t x t t x t x t x t t++++++=+∆=+∆=+∆ （2.20）2.3 微分方程的数值解法在分子动力学中的应用分子动力学在计算过程中需要计算每个粒子的运动状态，因此需要一套完整的迭代法则来计算每个粒子在每个时刻的位移、速度和加速度。

现在假设只有一个粒子，只考虑这个粒子在重力场内的自由落体运动状态。

因此问题描述为：()()()21()0 9.8/,1,1,40 2.0 00.0y t y t g e g m s m kg J m my y βαβαβ--+-======= （2.21）以下我们将利用不同的迭代方法解决上述问题并对比各种迭代方法进行对比。

2.3.1 欧拉法计算实例利用降阶的方法，将原问题（3.1）式化简为一个一次微分方程组。

设()()xt y t =因此可以得到：()()()()y t x t y t x t eg m βαβ-⎧=⎪⎨=-⎪⎩（2.22） 对上式两个一次微分项分别用差值逼近得到两个差分格式()()()()()()()()y t y t t y t x t y t tx t t x t x t e g m t βαβ-+∆-⎧==⎪⎪∆⎨+∆-⎪=-=⎪∆⎩（2.23） 将（3.3）式带入到（3.2）式就得到了原问题的欧拉法迭代格式()()()()()()y t x t t t e g x t m y t t ty t y t βαβ-⎧⎛⎫+∆=∆-+⎪⎪⎝⎭⎨⎪+∆=∆+⎩（2.24） 将上述欧拉法迭代格式用Fortran 程序实现为图2.1 欧拉法Fortran 程序利用欧拉法得到速度和位移曲线为v e l oc i t y （m /s ）time （s ）D i s p l a c e m e n t （m ）time （s ）图2.2 欧拉法速度-时间曲线 图2.3 欧拉法位移-时间曲线2.3.2 Newmark-beta 算法计算实例由于Newmark-beta 算法是利用加权的思想对二次微分项即加速度进行了修正，因此我们对原问题做出相应的变换，将原问题改写为力学运动方程常见的形式：()()()21() 9.8/,1,1,40 2.0 00.0y t y t e g g m s m kg J m my y βαβαβ--=-====== （2.25）将上式带入到Newmark-beta 算法的迭代格式（2.13）中去，并用Fortran 程序实现为图2.4 Newmark-beta 算法Fortran 程序 利用Newmark-beta 算法得到速度和位移曲线为v e l o c i ty (m /s )time(s)d i s p l a ce m e n t (m )time(s)图2.5 Newmark-beta 算法速度-时间曲线 图2.6 Newmark-beta 算法位移-时间曲线2.3.3 Leap-Frog （蛙跳）算法就算实例Leap-Frog （蛙跳）算法是将一次微分项即速度项进行了修正，它巧妙的利用二分之一迭代步的速度。