《数值计算：微分方程求解》由会员分享，可在线阅读，更多相关《数值计算：微分方程求解（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)
12、,(, ),(11cpnpnncnnnpyyyyxhfyyyxhfyy 解解：改进的：改进的Euler公式为公式为 )()()(cpnpnpncnnnnpyyyyxyhyyyxyhyy212211例例 作为比较，我们仍用作为比较，我们仍用Euler法中的那个例子。法中的那个例子。 1)0(2yyxyy)10( x1.73211.73791.01.67331.67820.91.61251.61530.81.54921.55250.71.48321.48600.61.41421.41640.51.34161.34340.41.26491.26620.31.18321.18410.21.09541.
13、09590.1nxny)(nxy我们仍取步长我们仍取步长h=0.1=0.1，计算结果，计算结果见右表：见右表：1 Eulers MethodMatlab作图显示作图显示 从表和图可以看出，改进的从表和图可以看出，改进的Euler法的精度提高了法的精度提高了不少。不少。1 Eulers Method初值问题的单步法可用一般形式表示为初值问题的单步法可用一般形式表示为: : 局部截断误差和方法的阶局部截断误差和方法的阶),(11hyyxhyynnnnn 称为增量函数称为增量函数 ),(),(111 nnnnnyxfhyyx 隐式欧拉法有隐式欧拉法有例如例如对欧拉法有对欧拉法有 ， ),(),(nn
14、nnyxfhyx 含有含有 时，方法是时，方法是隐式隐式的，的， 不含有不含有 时，方法是时，方法是显式显式的。的。1 ny 1 ny 从从 开始计算，如果考虑每一步产生的误差，直到开始计算，如果考虑每一步产生的误差，直到 则有误差则有误差 ，称为该方法在，称为该方法在 的的整体截断误整体截断误差，差，分析和求得整体截断误差是复杂的。为此，我们仅考分析和求得整体截断误差是复杂的。为此，我们仅考虑从虑从 到到 的局部情况，并假设的局部情况，并假设 之前的计算没有误之前的计算没有误差，即差，即 , ,下面给出单步法的下面给出单步法的局部截断局部截断误差概念。误差概念。nnnyxye )(0 xnx
15、nxnx1 nxnx)(nnxyy 1 Eulers Method即即 在假设在假设 yn = y(xn)，即第，即第n 步计算是精确的前提下，考步计算是精确的前提下，考虑的截断误差虑的截断误差 Tn+1= y(xn+1) yn+1 称为称为局部截断误差局部截断误差 /* local truncation error */。定义定义若某算法的局部截断误差为若某算法的局部截断误差为O(hp+1)，则称该算法有，则称该算法有p 阶精度。阶精度。 欧拉法的局部截断误差：欧拉法的局部截断误差：),()()()()()(321112nnnnhnnnnnyxhfyhOxyxyhxyyxyT )()(322
16、hOxynh 欧拉法具有欧拉法具有 1 阶精度。阶精度。Tn+1 的的主项主项/* leading term */为显式单步法的为显式单步法的局部截断误差局部截断误差。),(,()()(11hxyxhxyxyTnnnnn 定义定义 设设 是初值问题的准确解，称是初值问题的准确解，称)(xy即两步欧拉公式具有即两步欧拉公式具有 2 阶精度。阶精度。)()(3111hOyxyTnnn 两步两步欧拉法的局部截断误差：欧拉法的局部截断误差：1 Eulers Method11111)()4)(3)2)()1 nnnnnnnnnyxyTxxyxyxyy比比较较处处展展开开；在在将将准准确确解解）处处展展开
17、开；在在将将差差分分解解；作作局局部部化化假假设设，即即设设：求求局局部部截截断断误误差差的的步步骤骤)()(2)()()()()(2)()(,()()(322321111hxyhhxyhxyhhxyhxyhxyxhfxyxyTnnnnnnnnnn 例例 求隐式欧拉格式求隐式欧拉格式 的的局部截断误差。局部截断误差。),(111 nnnnyxhfyyTn+1 的的主项主项/* leading term */具有具有 1 阶精度。阶精度。例例 求梯形格式求梯形格式 的的局部局部截断误差。截断误差。),(),(2111 nnnnnnyxfyxfhyy以上定义对隐式单步法也是适用的。以上定义对隐式单
18、步法也是适用的。)()(12)()(2)()()(2)(!3)(2)()()(2)()(434232111hxyhhxyhxyhxyxyhxyhxyhxyhxyxyhxyxyTnnnnnnnnnnnnn 具有具有 2 阶精度。阶精度。Tn+1 的的主项主项/* leading term */ Hey! Isnt the leading term of the local truncation error of Eulers method ? Seems that we can make a good use of it )(22ihxy 1 Eulers Method方方 法法 1 Euler
19、s Method显式欧拉显式欧拉隐式欧拉隐式欧拉梯形公式梯形公式两步欧拉公式两步欧拉公式简单简单精度低精度低稳定性最好稳定性最好精度低精度低, 计算量大计算量大精度提高精度提高计算量大计算量大精度提高精度提高, 显式显式多一个初值多一个初值, 可能影响精度可能影响精度 Cant you give me a formula with all the advantages yet without any of the disadvantages? Do you think it possible? Well, call me greedy OK, lets make it possible.显然，
20、显然， , , 是在点是在点 , , 处的斜率。以上公式用到处的斜率。以上公式用到的是两个点的斜率的加权平均，它为构造算法提供了新的是两个点的斜率的加权平均，它为构造算法提供了新的途径。的途径。Runge-Kutta方法就是这种思想的体现和发展。方法就是这种思想的体现和发展。1k2knx1nx2 龙格龙格 - 库塔法库塔法 /* Runge-Kutta Method */建立高精度的单步递推格式。建立高精度的单步递推格式。单步递推法的单步递推法的基本思想基本思想是从是从 ( xn , yn ) 点出发，以点出发，以某一某一斜率斜率沿直线达到沿直线达到 ( xn+1 , yn+1 ) 点。欧拉法
21、及其各种变形点。欧拉法及其各种变形所能达到的最高精度为所能达到的最高精度为2阶阶。 考察改进的欧拉法，可以将其改写为：考察改进的欧拉法，可以将其改写为：),(),(2121121211hKyhxfKyxfKKKhyynnnnnn 斜率斜率一定取一定取K1、 K2 的的平均值平均值吗？吗？步长一定是一个步长一定是一个h 吗？吗？2 Runge-Kutta Method)(,()()(),(),()()(111 yhfxyxyxxyhxyxynnnnnn 2 Runge-Kutta Method首先希望能确定系数首先希望能确定系数 1、 2、p，使得到的算法格式有，使得到的算法格式有2阶阶精度，即
22、在精度，即在 的前提假设下，使得的前提假设下，使得 )()(3111hOyxyTnnn )(nnxyy Step 1: 将将 K2 在在 ( xn , yn ) 点作点作 Taylor 展开展开)(),(),(),(),(2112hOyxfphKyxphfyxfphKyphxfKnnynnxnnnn )()()(2hOxyphxynn Step 2: 将将 K2 代入第代入第1式，得到式，得到 )()()()()()()()(322212211hOxyphxyhyhOxyphxyxyhyynnnnnnnn 将改进欧拉法推广为：将改进欧拉法推广为：),(),(12122111phKyphxfKy
23、xfKKKhyynnnnnn 121 10 p),(),(),(),(),(),()(yxfyxfyxfdxdyyxfyxfyxfdxdxyyxyx 2 Runge-Kutta MethodStep 3: 将将 yn+1 与与 y( xn+1 ) 在在 xn 点的点的泰勒泰勒展开作比较展开作比较)()()()(322211hOxyphxyhyynnnn )()(2)()()(321hOxyhxyhxyxynnnn 要求要求 ，则必须有：，则必须有：)()(3111hOyxyTnnn 21,1221 p 这里有这里有 个未知个未知数，数， 个方程。个方程。32存在存在无穷多个解无穷多个解。所有满
24、足上式的格式统称为。所有满足上式的格式统称为2阶龙格阶龙格 - 库库塔格式塔格式。21, 121 p注意到，注意到， 就是改进的欧拉法。就是改进的欧拉法。 Q: 为获得更高的精度，应该如何进一步推广？为获得更高的精度，应该如何进一步推广？),(1nnyxfk . Liiinnkhyy11 )(,(232131333kakahcyhcxfknn 其中其中 , , 均为待定系数，确定这些系数均为待定系数，确定这些系数的步骤与前面相似。的步骤与前面相似。 11 Lii 111ijija1 ic),(11 ijjijininikahcyhcxfkLi, 3 , 2 ),(1222hkcyhcxfknn
25、 Runge-Kutta方法的一般形式方法的一般形式2 Runge-Kutta Method2 Runge-Kutta Method 最常用为四级最常用为四级4阶阶经典龙格经典龙格-库塔法库塔法 /* Classical Runge-Kutta Method */ ：),(),(),(),()22(34222312221432161hKyhxfKKyxfKKyxfKyxfKKKKKyynnhnhnhnhnnnhnn 1)0(2yyxyy) 10( x我们仍用前面的例子来看看四阶我们仍用前面的例子来看看四阶Runge-Kutta方法的效果。方法的效果。解解: : 对于本题，经典的四对于本题，经典
26、的四阶阶Runge-Kutta方法具有以下形式：方法具有以下形式： 334223112143211)(22222222)22(6hkyhxhkykkhyhxkhykkhyhxkhykyxykkkkkhyynnnnnnnnnnnnnn2 Runge-Kutta Method1.73211.73211.01.61251.61250.81.48321.48330.61.34161.34170.41.18321.18320.2nxny)(nxy 这里，我们取步长这里，我们取步长h= =0.2，下面是计算结果，下面是计算结果（符号的意义同前）（符号的意义同前） 从上图可以看出，每一个数值解都从上图可以看
27、出，每一个数值解都“准确准确”的落的落在了在了真实解真实解的曲线上。与改进的的曲线上。与改进的Euler法所做出来的图法所做出来的图比较，似乎精确性没有很明显的提高，但是不要忘了比较，似乎精确性没有很明显的提高，但是不要忘了在在用用Runge-Kutta方法时，我们取的步长是方法时，我们取的步长是h=0.2=0.2。实。实际上，际上，Runge-Kutta方法的精确性是方法的精确性是要高很多要高很多。2 Runge-Kutta MethodMatlab作图2 Runge-Kutta Method注：注： 龙格龙格-库塔法库塔法的主要运算在于计算的主要运算在于计算 Ki 的值，即计算的值，即计算
28、 f 的的值。值。Butcher 于于1965年给出了计算量与可达到的最高精年给出了计算量与可达到的最高精度阶数的关系：度阶数的关系：753可达到的最高精度可达到的最高精度642每步须算每步须算Ki 的个数的个数)(2hO)(3hO)(4hO)(5hO)(6hO)(4hO)(2nhO8n 由于龙格由于龙格-库塔法的导出基于泰勒展开，故精度主要受库塔法的导出基于泰勒展开，故精度主要受解函数的光滑性影响。对于光滑性不太好的解，最好解函数的光滑性影响。对于光滑性不太好的解，最好采用采用低阶算法低阶算法而将步长而将步长h 取小取小。HW: p.152-153 #2, # 3, # 43 收敛性与稳定性
29、收敛性与稳定性 /* Convergency and Stability */ 收敛性收敛性 /* Convergency */定义定义 若某算法对于任意固定的若某算法对于任意固定的 x = xn= x0 + n h，当，当 h0 ( 同时同时 n ) 时有时有 yn y( xn )，则称该算法是，则称该算法是收敛收敛的。的。 例：例：就初值问题就初值问题 考察欧拉显式格式的收敛性。考察欧拉显式格式的收敛性。 0)0(yyyy 解：解：该问题的精确解为该问题的精确解为 xeyxy 0)( 欧拉公式为欧拉公式为nnnnyhyhyy)1 (1 0)1 (yhynn 对任意固定的对任意固定的 x =
30、 xn = n h ，有，有nnxhhxnhyhyy )1()1(/10/0 )(0nxxyeyn ehhh /10)1(lim3 Convergency and Stability注：注：1、判断显式单步格式的收敛性，归结为验证、判断显式单步格式的收敛性，归结为验证增量函数增量函数 是否满足是否满足Lipschitz 条件。条件。 2、单步格式的整体截断误差由初值误差及局部截断误差决定，单步格式的整体截断误差由初值误差及局部截断误差决定，整整体截断误差体截断误差比比局部截断误差的阶数局部截断误差的阶数低一阶低一阶。 3、要构造、要构造高精度高精度的计算方法，只需设法的计算方法，只需设法提高局
31、部截断误差的阶提高局部截断误差的阶即即可。可。 （2）微分方程的初值是精确的）微分方程的初值是精确的. 。)()(pnnnhyxye 对于一个对于一个p 阶的显式单步法，若满足如下条件阶的显式单步法，若满足如下条件定理定理（1）增量函数）增量函数 关于关于y 满足满足Lipschitz条件，即存在常数条件，即存在常数 0 L，使，使 RyyyyLhyxhyx ,),(),( 成立；成立；3 Convergency and Stability更详细的请参看教材更详细的请参看教材p137则该方法收敛，其整体截断误差为则该方法收敛，其整体截断误差为:3 Convergency and Stabili
32、ty 稳定性稳定性 /* Stability */例：例：考察初值问题考察初值问题 在区间在区间0, 0.5上的解。上的解。分别用欧拉显、隐式格式和改进的欧拉格式计算数值解。分别用欧拉显、隐式格式和改进的欧拉格式计算数值解。 1)0()(30)(yxyxy0.00.10.20.30.40.5精确解精确解改进欧拉法改进欧拉法 欧欧拉隐式拉隐式欧拉显式欧拉显式 节点节点 xixey30 1.0000 2.0000 4.0000 8.0000 1.6000 101 3.2000 101 1.00002.5000 10 1 6.2500 10 21.5625 10 23.9063 10 39.7656
33、 10 41.00002.50006.25001.5626 1013.9063 1019.7656 1011.00004.9787 10 22.4788 10 31.2341 10 46.1442 10 63.0590 10 7What is wrong ?! An Engineer complains: Math theorems are so unstable that a small perturbation on the conditions will cause a crash on the conclusions!3 Convergency and Stability定义定义若某算
34、法在计算过程中任一步产生的误差在以后的计若某算法在计算过程中任一步产生的误差在以后的计算中都算中都逐步衰减逐步衰减，则称该算法是，则称该算法是绝对稳定的绝对稳定的 /*absolutely stable */。一般分析时为简单起见，只考虑一般分析时为简单起见，只考虑试验方程试验方程 /* test equation */0Re( ) yy常数，可以常数，可以是复数是复数当步长取为当步长取为 h 时，将某算法应用于上式，并假设只在初值时，将某算法应用于上式，并假设只在初值产生误差产生误差 ，则若此误差以后逐步衰减，就称该，则若此误差以后逐步衰减，就称该算法相对于算法相对于 绝对稳定绝对稳定， 的
35、全体构成的全体构成绝对稳定区域绝对稳定区域。我们称我们称算法算法A 比算法比算法B 稳定稳定，就是指，就是指 A 的绝对稳定区域比的绝对稳定区域比 B 的的大大。000yy h h h3 Convergency and Stability例：例：考察显式欧拉法考察显式欧拉法011)1(yhyhyyiiii 000yy 011)1(yhyii 01111)1( iiiihyy由此可见，要保证初始误差由此可见，要保证初始误差 0 以后逐步衰减，以后逐步衰减，必须满足：必须满足：hh 1|1| h0-1-2ReIm(h)例：例：考察隐式欧拉法考察隐式欧拉法11 iiiyhyy iiyhy 11101
36、111 iih可见绝对稳定区域为：可见绝对稳定区域为：1|1| h210ReIm(h)注：注：一般来说，隐式欧拉法的绝对稳定性比同阶的显式一般来说，隐式欧拉法的绝对稳定性比同阶的显式法的好。法的好。3 Convergency and Stability例：例：隐式龙格隐式龙格-库塔法库塔法 ),., 1().,(.11111mjhKhKyhxfKKKhyymmjjijijmmii 而而显式显式 1 4 阶方法的绝对稳定阶方法的绝对稳定区域为区域为 )2,2(1111KhyhxfKhKyyiiii其中其中2阶方法阶方法 的绝对稳定区域为的绝对稳定区域为0ReIm(h)k=1k=2k=3k=4-1-2-3-123ReIm(h)无条件稳定无条件稳定HW: p.153 #6,8
