数值计算：微分方程求解

《数值计算：微分方程求解》由会员分享，可在线阅读，更多相关《数值计算：微分方程求解（29页珍藏版）》请在文档大全上搜索。

1、第五章第五章 常微分方程数值解常微分方程数值解/* Numerical Methods for Ordinary Differential Equations */ 考虑考虑一阶一阶常微分方程的常微分方程的初值问题初值问题 /* Initial-Value Problem */: 0)(,),(yaybaxyxfdxdy只要只要 f (x, y) 在在a, b R1 上连续，且关于上连续，且关于 y 满足满足 Lipschitz 条件条件，即存在与即存在与 x, y 无关的常数无关的常数 L 使使对任意定义在对任意定义在 a, b 上的上的 y1(x) 和和 y2(x) 都成立，则上述都成立，

2、则上述IVP存存在唯一解在唯一解。| ),(),(|2121yyLyxfyxf 本章的任务：计算出解函数本章的任务：计算出解函数 y(x) 在一系列节点在一系列节点 a = x0 x1 xn= b 处的近似值处的近似值 。 ),., 1()(Nnxyynn 建立常微分方程数值方法的基本思想建立常微分方程数值方法的基本思想 微分方程数值解法，其实是求出方程的解微分方程数值解法，其实是求出方程的解 在一系在一系列离散点上的近似值。则微分方程数值解的基本思想是：列离散点上的近似值。则微分方程数值解的基本思想是：求解区间和方程离散化。求解区间和方程离散化。)(xy 将求解区间将求解区间 离散化离散化,

3、 ,是在是在 上插入一系列的分上插入一系列的分点点 , ,使使bxxxxxaNnn .110 ba, ba, kx，求解区间离散化求解区间离散化 记记 称为步长，一般取称为步长，一般取hn = h 称为等步长节点称为等步长节点 。)1,.,1 , 0(1 Nnxxhnnnnhxxn 0),.,2, 1 ,0(Nn Nabh (常数常数),节点为节点为 mnxxnmdxxyxfxyxy)(,()()()(00yxy 然后将上式右端采用第四章介绍的数值积分离散化，然后将上式右端采用第四章介绍的数值积分离散化，从而获得原初值问题的一个离散差分格式。从而获得原初值问题的一个离散差分格式。将微分方程离散

4、化将微分方程离散化 将微分方程离散化，通常有下述方法：将微分方程离散化，通常有下述方法：（1 1）差商逼近法差商逼近法 即是用适当的差商逼近导数值。即是用适当的差商逼近导数值。（2 2）数值积分法数值积分法 基本思想是先将问题转化为积分方程基本思想是先将问题转化为积分方程（3）Taylor展开法展开法 见后面的叙述。见后面的叙述。1 欧拉方法欧拉方法 /* Eulers Method */ 欧拉公式：欧拉公式：x0 x1向前差商近似导数向前差商近似导数hxyxyxy)()()(010 ),()()()(000001yxfhyxyhxyxy 1y记为记为)1,., 0(),(1 Nnyxfhyy

5、nnnn亦称为亦称为欧拉折线法欧拉折线法 /* Eulers polygonal arc method*/ 1)0(2yyxyy)10( x例例. . 求初值问题求初值问题解解：本题的：本题的Euler公式的具体形式为公式的具体形式为)(nnnnnyxyhyy21 1 Eulers Method1.73211.78481.01.41421.43510.51.67331.71780.91.24161.35820.41.61251.64980.81.26491.27740.31.54921.58030.71.18321.19180.21.48321.50900.61.09541.10000.1nx

6、nxnyny)(nxy)(nxy 如果说表格仍不够直观的话，我们用如果说表格仍不够直观的话，我们用MatlabMatlab做出积分曲做出积分曲线与近似值的图，如下：线与近似值的图，如下：取步长取步长h=0.1=0.1。我们将计算结果与其解析解的精确值一同。我们将计算结果与其解析解的精确值一同列在下表中，其中列在下表中，其中 是节点，是节点， 是节点上的近似值，是节点上的近似值，是精确值，结果见下表：是精确值，结果见下表：nynx)(nxy1 Eulers Method 从图中可以看出，灰色连续的曲线就是初值问题的解析解从图中可以看出，灰色连续的曲线就是初值问题的解析解 的曲线。而红色的星点便是

7、数值解。在图上似乎的曲线。而红色的星点便是数值解。在图上似乎数值解与曲线的偏差不是很大，但不要忘记这只是在数值解与曲线的偏差不是很大，但不要忘记这只是在0 0到到1 1范围范围内的。通过后面用其他方法解本题，大家便会发现内的。通过后面用其他方法解本题，大家便会发现EulerEuler方法误方法误差其实是很大的。差其实是很大的。xy21MatlabMatlab作图显示作图显示1 Eulers Method1 Eulers Method 隐式欧拉法隐式欧拉法 /* implicit Euler method */向后差商近似导数向后差商近似导数hxyxyxy)()()(011 x0 x1)(,()

8、(1101xyxfhyxy ),( y1( fh11.0,1 Nnxn yynnn) 由于未知数由于未知数 yn+1 同时出现在等式的两边，同时出现在等式的两边，不能直接得到，故称为不能直接得到，故称为隐式隐式 /* implicit */ 欧拉公式，而前者称欧拉公式，而前者称为为显式显式 /* explicit */ 欧拉公式。欧拉公式。一般先用显式计算一个初值，再一般先用显式计算一个初值，再迭代迭代求解。即求解。即),()0(1nnnnyxhfyy ),()(11)1(1knnnknyxhfyy ,2,1 ,0 k如果迭代过程收敛，则某步后如果迭代过程收敛，则某步后 就可以作为就可以作为

9、，从而进行，从而进行下一步的计算。下一步的计算。)(1kny 1 ny 梯形方法的平均化思想梯形方法的平均化思想可以借助几何直观说明，同可以借助几何直观说明，同Euler方法的图，见右：方法的图，见右：1 Eulers Method 梯形公式梯形公式 /* trapezoid formula */ 显、隐式两种算法的显、隐式两种算法的平均平均)1,., 0(),(),(2111 Nnyxfyxfhyynnnnnn 两步欧拉公式两步欧拉公式 /* midpoint formula */中心差商近似导数中心差商近似导数hxyxyxy2)()()(021 x0 x2x1)(,(2)()(1102xy

10、xfhxyxy 1,., 1),(211 NnyxfhyynnnnnPABnx1nxQ需要需要2个初值个初值 y0和和 y1来启动递推来启动递推过程，这样的算法称为过程，这样的算法称为两步法两步法 /* double-step method */，而前面的三种算法都是，而前面的三种算法都是单步法单步法 /* single-step method */。1 nP 改进欧拉法改进欧拉法 /* modified Eulers method */Step 1: 先用先用显式显式欧拉公式作欧拉公式作预测预测，算出，算出),(1nnnnyxfhyy Step 2: 再将再将 代入代入隐式隐式梯形公式的右边

11、作梯形公式的右边作校正校正，得到，得到1 ny),(),(2111 nnnnnnyxfyxfhyy注：注：此法亦称为此法亦称为预测预测-校正法校正法 /* predictor-corrector method */。可以证明该算法具有可以证明该算法具有 2 阶精度，同时可以看到它是个阶精度，同时可以看到它是个单单步步递推格式，比隐式公式的迭代求解过程递推格式，比隐式公式的迭代求解过程简单简单。后面将。后面将看到，它的看到，它的稳定性高稳定性高于显式欧拉法。于显式欧拉法。 )1,., 0(),(,),(211 Nnyxfhyxfyxfhyynnnnnnnn1 Eulers Method)(21)