偏微分方程的几种数值解法及其应用

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

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

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

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

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

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

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

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

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

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

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

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

()
()00(,)dx t f x t dt
x 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 t
dx t x t x t dt t
dx t x t x t dt t
++++--=∆-=
∆-=∆ （2.2） 上式中，分别为向前差分法、向后差分法、中心差分法。

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

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

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

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

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

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

根据运动学方程：
2
12
x 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）
同理，对原函数也进行相同的处理，利用广义中值定理对位移进行修正得到：
()()()()211
2
l 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 算法的迭代格式为
()()()()()()()()()()()()()
1112
11
2
2
2
l l l l l l l l l l x t x t x t t
x 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/22
l l l x t x t x t -++=
（2.16）
在迭代启动的时候需要对初始条件进行修改，需要用初始条件来反算()1/2x
t -：
()()()1/200/2x t x t x t t -=-∆ （2.17）
2.2.5 Velocity Verlet 算法
Velocity Verlet 算法由于其速度计算精度更高被广泛用在分子动力学计算中，它的基本原理为，对位移进Taylor 展开，其展开式为：
()()()()()23111
26
l l l l i l x t x t x t t x t t b t t +=+∆+
∆+∆ （2.18） 略去高阶项就得到了Velocity Verlet 算法的迭代格式
()()()()211
2
l l l l x t x t x t t x t t +=+∆+∆ （2.19a ） ()()()()111
2l l l l x t x t x t x t t ++=+
+∆⎡⎤⎣
⎦ （2.19b ） 还等价于
()()()()()()()()()1/211/211/211
2
1
2
l l l l l l l l l x t x t x t t
x 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.0
y t y t g e g m s m kg J m m
y y βαβ
αβ--+-
======= （2.21）
以下我们将利用不同的迭代方法解决上述问题并对比各种迭代方法进行对比。

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

设()()x
t y t =因此可以得到：
()()()()
y t x t y t x t e
g m βαβ-⎧=⎪⎨=-⎪⎩
（2.22） 对上式两个一次微分项分别用差值逼近得到两个差分格式
()()()()()()()()y t y t t y t x t y t t
x 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 o
c 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.0
y t y t e g g m s m kg J m m
y y βαβ
αβ--=
-====== （2.25）
将上式带入到Newmark-beta 算法的迭代格式（2.13）中去，并用Fortran 程序实现为
图2.4 Newmark-beta 算法Fortran 程序 利用Newmark-beta 算法得到速度和位移曲线为
v e l o c i t
y (m /s )
time(s)
d i s p l a c
e m e n t (m )
time(s)
图2.5 Newmark-beta 算法速度-时间曲线 图2.6 Newmark-beta 算法位移-时间曲线
2.3.3 Leap-Frog （蛙跳）算法就算实例
Leap-Frog （蛙跳）算法是将一次微分项即速度项进行了修正，它巧妙的利用二分之一迭代步的速度。

因此将修改后的原问题公式（3.5）式带入到Leap-Frog （蛙跳）算法迭代格式中，就得到了原问题的Leap-Frog （蛙跳）算法解答。

并利用Fortran 程序实现其迭代过程如下：
图2.7 Leap-Frog （蛙跳）算法Fortran 程序 利用Leap-Frog （蛙跳）算法得到速度和位移曲线为
v e l o c i t
y (m /s )
time(s)
d i s p l a c
e m e n t (m )
time(s)
图2.8 Leap-Frog （蛙跳）算法速度-时间曲线 图2.9 Leap-Frog （蛙跳）算法位移-时间曲线
2.3.4 Velocity Verlet 算法计算实例
Velocity Verlet 算法对速度的计算结果较精确，被广泛应用在分子动力学计算中，其核心思想是将位移函数进行Talor 函数展开，因此将修改后的原问题带入到Velocity Verlet 算法的迭代格式中，就得到原问题的Velocity Verlet 算法解答。

并利用Fortran 程序实现其迭代过程如下：
图2.10 Velocity Verlet 算法Fortran 程序 利用Velocity Verlet 算法得到的速度和位移曲线为
v e l o c i t
y (m /s )
time(s)
d i s p l a c
e m e n t (m )
time(s)
图2.11 Velocity Verlet 算法速度-时间曲线 图2.12 Velocity Verlet 算法位移-时间曲线
2.3.5 其他工况下的欧拉迭代方法计算
利用欧拉迭代的方法对上抛运动、有阻尼运动这两种种工况分别做了计算，计算结果和实际物理过程相符。

1）上抛运动
上抛运动和自由落体运动相比，差别在于其初始速度的差别，自由落体的初始速度为零，而上抛运动的初始速度为一个向上的速度。

因此只修改原问题的初始速度就得到上抛运动的问题描述。

()()()()21
()20 9.8/,1,1,40 2.0 00.5
y t y t g e y t g m s m kg J m m
y y βαβ
αβ--+-
+======= （2.26）
将上述问题带入到欧拉法的迭代格式中，得到速度和位移曲线。

v e l o c i t y
(m /s )
time （s ）
d i s p l a c
e m e n t (m )
time(s)
图2.11 上抛运动速度-时间曲线 图2.12 上抛运动位移-时间曲线 2） 有阻尼运动
有阻尼运动和自由落体运动相比，差别在于有阻尼运动在原自由落体的问题中多产生一项阻尼项，修改原问题，在原问题的基础上添加一项阻尼项得到有阻尼运动的问题描述。

()()()()21
()20 9.8/,1,1,40 2.0 00.0
y t y t g e y t g m s m kg J m m
y y βαβ
αβ--+-
+======= （2.27）
将上述问题带入到欧拉迭代格式中，得到速度和位移曲线。

v e l o c i t y (m /s )
time(s)
d i s p l a c
e m e n t (m )
time(s)
图2.11 有阻尼运动速度-时间曲线 图2.12 有阻尼运动位移-时间曲线
2.4 微分方程四种数值解法的对比
上文介绍了微分方程的四种数值解法欧拉法，Newmark-beta 算法、Leap-Frog （蛙跳）算法和Velocity Verlet 算法。

下面我们针对自由落体运动工况和有阻尼运动工况对四种数值解法进行对比。

在自由落体的情况下，分别用四种数值迭代方法求解自由运动方程得到速度和位移的曲线。

v
e
l
o
c
i
t
y
(
m
/
s
)
time(s)
d
i
s
p
l
a
c
e
m
e
n
t
(
m
)
time(s)
图2.13 四种数值解法速度对比图图2.14 四种数值解法位移对比图通过自由落体工况下，四种数值解法的对比可以得出，四种数值解法在求解自由落体这种工况下速度计算结果差距不大，速度计算结果从大到小排列为Euler>Leap-Frog>Newmark=Velocity Verlet，位移结果上差距也不是很明显，位移计算结果从大到小排列为Euler>Newmark>Velocity Verlet>Leap-Frog。

在有阻尼运动的工况下，分别用四种数值迭代方法来求解有阻尼运动方程得到速度和位移的曲线。

v
e
l
o
c
i
t
y
(
m
/
s
)
time(s)
d
i
s
p
l
a
c
e
m
e
n
t
(
m
)
time(s)
图2.13 四种数值解法速度对比图图2.14 四种数值解法位移对比图在有阻尼运动的工况下，对四种数值方法对比后发现，Leap-Frog（蛙跳）算法的收敛速度较快，大约在1750步左右就收敛，其原因是Leap-Frog（蛙跳）算法是对一次微分项进行修正，Leap-Frog（蛙跳）算法每一个迭代步的速度是由两个二分之一迭代步速度的平均得到的，因此在后一个二分之一迭代步的速度逐渐收敛后，速度项会快速收敛。

而其他算法都是对二次微分项进行修正，当二次微分项逐渐收敛时，速度项才开始逐渐收敛。

这就是Leap-Frog（蛙跳）算法收敛速度比其他算法收敛快的根本原因。

从对比过程中我们们发现欧拉法虽然原理简单使用起来比较方便，但是欧拉法算的速度项偏大，而且收敛速度最慢，因此不适合做大规模迭代；Leap-Frog（蛙跳）算法收敛速度过慢，速度计算不是很精确，Newmark-beta算法一般在有限元动力学分析中应用广泛。

Velocity Verlet算法收敛速度适中，对速度的计算较精确，因此Velocity Verlet算法在分子动力学计算中被广泛的应用。

3 总结
本文总结了四种常用微分方程的数值解法，并较深入的剖析了四种数值解法的基本原理。

运用微分方程的四种不同解法对一个粒子在重力场发生碰撞反弹过程进行计算并对比四种不同算法得到的结果。

经过对比后发现，Newmark-beta算法和Velocity Verlet算法在收敛速度和结果的精确性上占有较大优势，欧拉法收敛速度过慢，L eap-Frog（蛙跳）算法收敛速度过快，对速度的计算结果精度过于粗糙。

Newmark-beta 算法在有限元动力分析中应用广泛，Velocity Verlet算法分子动力学计算中应用广泛。