《偏微分方程的几种数值解法及其应用》由会员分享，可在线阅读，更多相关《偏微分方程的几种数值解法及其应用（12页珍藏版）》请在金锄头文库上搜索。
1、1 1 1 常微分常微分方程及其方程及其数值数值解法解法 1 1.1 .1 常微分常微分方程方程概述概述 在数学上，物质的运动和变化规律是通过函数关系来表示的，在一些复杂的现象中，我们要求的未知 量就变成了满足特定条件的一个或一些未知函数。有的时候，我们需要利用导数或者微分的关系，即这些 未知函数的导数与自变量满足某种关系，这种方程我们称之为微分方程。未知函数是一元函数的微分方程 称之为常微分方程， 未知函数是多元函数的微分方程我们称之为偏微分方程， 我们这里只考虑常微分方程。 常微分方程的解，就是找出一个代入方程使之成为恒等式的函数。若微分方程的解中含有任意常数的 个数与方程阶数相同，且任意
2、常数之间不能合并，则称此解为该方程的通解。当通解中的各任意常数都取 特定值时所得到的解，称为方程的特解。在实际问题中，这些函数往往还需要满足一些特定条件，这称之 为定解条件。 但在实际问题中，很多常微分方程的解析表达式过于复杂，甚至得不到通解的解析表达式。而且，常 微分方程的特解是否存在，存在几个特解，这涉及到微分方程解的存在性和唯一性定理。因此，在实际应 用中， 我们通常利用数值的方法来求得方程的数值解， 在误差允许的范围内， 我们用数值解来替代解析解。 所以，研究常微分的数值解法是很有必要的。 2.2 2.2 常微分常微分方程的数值解法方程的数值解法 常微分方程的数值解法是有常微分方程的定
3、解条件提出的，首先我们考虑如下一阶常微分方程的初值 问题。 00 ( , ) dx t f x t dt x tx （2.1） 2.2.12.2.1 欧拉法欧拉法 欧拉法 （又称差分法） 是常微分方程初值问题数值解法中最简单最古老的方法， 它的基本思路是将 （2.1） 式中导数项用差分来逼近，从而将一个微分方程转化为一个代数方程，以便迭代求解。根据用于逼近的差 分方式来分，可以分为向前差分、向后差分、中心差分。 2 1 11 11 2 lll lll lll dx tx tx t dtt dx tx tx t dtt dx tx tx t dtt （2.2） 上式中，分别为向前差分法、向后差分
4、法、中心差分法。 以向前差分法为例，我们在这里推导向前欧拉法的迭代格式。由（2.2）式中第一个向前差分法公式我 们可以得出 1 l ll dx t x tx tt dt （2.3） 再由原问题（2.1）可以得到向前欧拉法的迭代格式为： 1 , lll x tx ttf x t （2.4） 因此我们可以利用初始条件，依次得到后一个时刻方程的解，直至求得的解收敛到所要求的时刻为止。 当原问题为更为复杂的常微分方程组或者高阶常微分时，只需要将x看做向量，原问题就可以简化或 者降阶为一个一阶常微分方程组或者一个一阶常微分方程。这样就可以利用欧拉法的迭代格式进行迭代计 算。 2.2.2 2.2.2 改进
5、改进欧拉法欧拉法 用数值积分的方法对原问题（2.1）式进行离散化处理，方程两边做积分有： 1 1 ( ), l l t ll t x tx tf x t t dt （2.5） 对右端积分使用梯形积分公式可得 1 11 ( ), 2 l l t llll t t f x t t dtf x ttf x tt （2.6） 将（2.6）式带入到（2.5）式中可以得到 111 , 2 llllll t x tx tf x ttf x tt （2.7） 在实际计算过程中，改进的欧拉法是用欧拉法先求得一个预测值 1l x t ，再利用这个预测值来来计算 1l x t ，因此改进欧拉法的迭代格式为： 3 1
6、 111 , , 2 llll llllll x tx ttf x tt t x tx tf x ttf x tt （2.8） 从（2.7）式不难看出改进欧拉法是隐式迭代算法，而欧拉法是显示迭代算法。 2.2.3 2.2.3 NewmarkNewmark- -betabeta 算法算法 Newmark-beta 算法是用于解决微分方程的数值积分算法之一， 它通常用于有限单元法中的动力学分析。 根据运动学方程： 2 1 2 xxtxt （2.7） 利用广义中值定理，Newmark-beta 算法对于一次微分项修改为 1ll x tx ttxt （2.8） 和欧拉法不同的是，Newmark-bet
7、a 算法对于二次微分项不是直接利用上一个迭代步中的二次微分结果，而 是对当前迭代步和上一迭代步两个迭代步的二次微分结果进行了加权取值，利用加权修正后的二次微分结 果进行迭代计算。 1 1( ) 01 ll xtx tx t （2.9） 将（2.9）式带入到（2.8）式中得到，一次微分项的迭代格式为： 11 1 llll x tx ttx ttx t （2.10） 同理，对原函数也进行相同的处理，利用广义中值定理对位移进行修正得到： 2 1 1 2 lll x tx ttx tt xt （2.11） 依旧利用加权来修正上式中的二次微分项为 1 1 22 01 ll xtx tx t （2.12）
8、 现假设 Newmark-beta 算法中令0,1/ 2得到，Newmark-beta 算法的迭代格式为 1 11 2 1 1 2 2 2 ll llll llll xtx tx t t x tx tx tx t t x tx ttx tx t （2.13） 由上式可以利用初始条件依次得到每个迭代步的数值解。 4 2 2.2.4 .2.4 LeapLeap- -FrogFrog（蛙跳）算法（蛙跳）算法 Leap-Frog（蛙跳）算法是对一次微分项即速度进行修正，它的基本原理是，首先利用当前时刻的加速 度，计算半个时间步长后的速度： 1/21/2lll x tx tx tt （2.14） 之后，
9、计算下一个步长时刻的位置： 11/2lll x tx tx tt （2.15） 利用两个半个时间步长的速度来计算当前时刻的速度 1/21/2 2 ll l x tx t x t （2.16） 在迭代启动的时候需要对初始条件进行修改，需要用初始条件来反算 1/2 x t： 1/200 /2x tx tx tt （2.17） 2.2.52.2.5 Velocity VerletVelocity Verlet 算法算法 Velocity Verlet 算法由于其速度计算精度更高被广泛用在分子动力学计算中，它的基本原理为，对 位移进 Taylor 展开，其展开式为： 23 1 11 26 llllil
10、 x tx tx ttx ttb tt （2.18） 略去高阶项就得到了 Velocity Verlet 算法的迭代格式 2 1 1 2 llll x tx tx ttx tt （2.19a） 11 1 2 llll x tx tx tx tt （2.19b） 还等价于 1/2 11/2 11/21 1 2 1 2 lll lll lll x tx tx tt x tx tx tt x tx tx tt （2.20） 2.32.3 微分方程的微分方程的数值解法在分子动力学中的应用数值解法在分子动力学中的应用 5 分子动力学在计算过程中需要计算每个粒子的运动状态，因此需要一套完整的迭代法则来计算
11、每个粒 子在每个时刻的位移、速度和加速度。现在假设只有一个粒子，只考虑这个粒子在重力场内的自由落体运 动状态。因此问题描述为： ( )21 ( )0 9.8/,1,1 ,4 02.0 00.0 y t y tgegm smkgJm m yy （2.21） 以下我们将利用不同的迭代方法解决上述问题并对比各种迭代方法进行对比。 2.3.2.3.1 1 欧拉法欧拉法计算计算实例实例 利用降阶的方法，将原问题（3.1）式化简为一个一次微分方程组。设 ( )x ty t因此可以得到： y t x ty t x teg m （2.22） 对上式两个一次微分项分别用差值逼近得到两个差分格式 y t y tt
12、y t x ty t t x ttx t x teg mt （2.23） 将（3.3）式带入到（3.2）式就得到了原问题的欧拉法迭代格式 y t x tttegx t m y ttty ty t （2.24） 将上述欧拉法迭代格式用 Fortran 程序实现为 6 图 2.1 欧拉法 Fortran 程序 利用欧拉法得到速度和位移曲线为 0.00.20.40.60.81.0 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 velocity（m/s） time（s） velocity 0.00.20.40.60.81.0 0.00 0.05 0.10 0.15 0
13、.20 Displacement（m） time（s） Displacement 图 2.2 欧拉法速度-时间曲线 图 2.3 欧拉法位移-时间曲线 2.3.2.3.2 2 NewmarkNewmark- -betabeta 算法算法计算实例计算实例 由于 Newmark-beta 算法是利用加权的思想对二次微分项即加速度进行了修正，因此我们对原问题做 出相应的变换，将原问题改写为力学运动方程常见的形式： ( )21 ( ) 9.8/,1,1 ,4 02.0 00.0 y t y teggm smkgJm m yy （2.25） 将上式带入到 Newmark-beta 算法的迭代格式（2.13）中去，并用 Fortran 程序实现为 7 图 2.4 Newmark-beta 算法 Fortran 程序 利用 Newmark-beta 算法得到速度和位移曲线为 0.00.10.20.30.40.5 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 velocity(m/s) ti
