《第四章常微分方程数值解法优秀课件.ppt》由会员分享，可在线阅读，更多相关《第四章常微分方程数值解法优秀课件.ppt（65页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、第四章常微分方程数值第四章常微分方程数值解法解法1 1第1页，本讲稿共65页常微分方程数值解法引言（常微分方程数值解法概述）引言（常微分方程数值解法概述）显式欧拉法、隐式欧拉法、二步欧拉法显式欧拉法、隐式欧拉法、二步欧拉法局部截断误差与精度局部截断误差与精度改进的欧拉方法改进的欧拉方法龙格龙格-库塔方法库塔方法收敛性与稳定性简述收敛性与稳定性简述一阶常微分方程组与高阶常微分方程一阶常微分方程组与高阶常微分方程2第2页，本讲稿共65页引言引言一阶常微分方程初值问题：一阶常微分方程初值问题：微分方程初始条件定理：若 f(x,y)在某闭区域 R：上上连连续续，且且在在 R 域域内内满满足足李李普普希
2、希兹兹(Lipschitz)条条件件，即即存存在在正正数数 L，使使得得对对于于 R 域域内内的的任任意意两两值值 y1,y2，下下列不等式成立：列不等式成立：则上述初值问题的连续可微的解则上述初值问题的连续可微的解 y(x)存在并且唯一。存在并且唯一。3第3页，本讲稿共65页引言（续）引言（续）实实际际生生产产与与科科研研中中，除除少少数数简简单单情情况况能能获获得得初初值值问问题题的的初初等等解解（用用初初等等函函数数表表示示的的解解）外外，绝绝大大多多数数情情况况下是求不出初等解的。下是求不出初等解的。有有些些初初值值问问题题即即便便有有初初等等解解，也也往往往往由由于于形形式式过过于于
3、复复杂而不便处理。杂而不便处理。实实用用的的方方法法是是在在计计算算机机上上进进行行数数值值求求解解：即即不不直直接接求求 y(x)的的显显式式解解，而而是是在在解解所所存存在在的的区区间间上上，求求得得一一系系列点列点 xn(n=0,1,2,)上解的近似值。上解的近似值。4第4页，本讲稿共65页欧拉欧拉（Euler）方法方法方法一方法一化导数为差商的方法化导数为差商的方法由于在逐步求解的过程中，由于在逐步求解的过程中，y(xn)的准确值无法求解的准确值无法求解出来，因此用其近似值代替。出来，因此用其近似值代替。为避免混淆，以下学习简记：为避免混淆，以下学习简记：y(xn)：待求函数：待求函数
4、 y(x)在在 xn 处的精确函数值处的精确函数值yn：待求函数：待求函数 y(x)在在 xn 处的近似函数值处的近似函数值5第5页，本讲稿共65页代入初值问题表达式可得：代入初值问题表达式可得：根据根据 y0 可以一步步计算出函数可以一步步计算出函数 y=y(x)在在 x1,x2,x3 x4,上的近似值上的近似值 y1,y2,y3,y4,常微分方程数值解是一组离散的函数值数据，它的常微分方程数值解是一组离散的函数值数据，它的精确表达式很难求解得到，但可以进行插值计算后精确表达式很难求解得到，但可以进行插值计算后用插值函数逼近用插值函数逼近 y(x)6第6页，本讲稿共65页欧拉方法（续）欧拉方
5、法（续）方法二方法二数值积分法数值积分法同样以近似值同样以近似值 yn 代替精确代替精确值值 y(xn)可得可得：将微分方程将微分方程 y =f(x,y)在区在区间间 xn,xn+1 上上积积分：分：7第7页，本讲稿共65页欧拉方法的几何意义欧拉方法的几何意义xy08第8页，本讲稿共65页隐式欧拉法隐式欧拉法在数值积分法推导中，积分的近似值取为积分区间宽在数值积分法推导中，积分的近似值取为积分区间宽度与右端点处的函数值乘积，即：度与右端点处的函数值乘积，即：这样便得到了隐式欧拉法：这样便得到了隐式欧拉法：含有未知的函数值隐式欧拉法没有显式欧拉法方便隐式欧拉法没有显式欧拉法方便9第9页，本讲稿共
6、65页二步欧拉法二步欧拉法在在数数值值积积分分法法推推导导中中，积积分分区区间间宽宽度度选选为为两两步步步步长长，即积分区间为：即积分区间为：xn-1,xn+1，则：，则：以以 y(x)在在 xn-1,xn 上的近似值上的近似值代替精确代替精确值值可得可得：需要前两步的计算结果中矩形公式中矩形公式10第10页，本讲稿共65页梯形公式欧拉法梯形公式欧拉法在数值积分法中，如果用梯形公式近似计算在数值积分法中，如果用梯形公式近似计算 f(x,y)在区间在区间 xn,xn+1 上的积分，即：上的积分，即：用近似值代替精确值可得梯形公式欧拉法：用近似值代替精确值可得梯形公式欧拉法：上式右端出现了未知项，
7、可见梯形法是隐式欧拉法上式右端出现了未知项，可见梯形法是隐式欧拉法的一种；实际上，梯形公式欧拉法是显式欧拉法与的一种；实际上，梯形公式欧拉法是显式欧拉法与隐式欧拉法的隐式欧拉法的算术平均算术平均。11第11页，本讲稿共65页例例用显式欧拉法、隐式欧拉法、梯形法求解初值问题：用显式欧拉法、隐式欧拉法、梯形法求解初值问题：取取 h=0.1，计算到，计算到 x=0.5，并与精确解进行比较，并与精确解进行比较解：由已知条件可得：解：由已知条件可得：h=0.1，x0=0，y0=1，f(x,y)=-=-y+x+1显式欧拉法：12第12页，本讲稿共65页例：（续）例：（续）隐式欧拉法：化简得：梯形公式欧拉法
8、：13第13页，本讲稿共65页计算结果：计算结果：xn显式法 yn隐式法 yn梯形法 yn精确解 y(xn)0.011110.11.0000001.0090911.0047621.0048370.21.0100001.0264461.0185941.0197310.31.0290001.0513151.0406331.0408180.41.0561001.0830141.0700971.0703200.51.0904901.1209221.1062781.106531本题的精确解为：显式法隐式法梯形法14第14页，本讲稿共65页局部截断误差局部截断误差为了简化分析某常微分方程数值算法的误差，现
9、假为了简化分析某常微分方程数值算法的误差，现假设设 yn=y(xn)，即，即在前一步在前一步 yn 准确的前提下准确的前提下，估计：，估计：称称上上述述误误差差 Tn+1 为为该该常常微微分分方方程程数数值值算算法法的的局局部部截截断误差断误差如如果果某某个个常常微微分分方方程程数数值值算算法法的的局局部部截截断断误误差差可可表表示为示为 O(h p+1)，则称该数值算法的精度是，则称该数值算法的精度是 p 阶阶欧欧拉拉法法的的精精度度为为一一阶阶；二二步步欧欧拉拉法法的的精精度度为为二二阶阶；梯形法的精度为二阶。梯形法的精度为二阶。15第15页，本讲稿共65页泰勒展开法泰勒展开法如果初值问题
10、中的如果初值问题中的 f(x,y)充分可微，则可将充分可微，则可将 y(xn+1)在点在点 xn 处展开：处展开：如果只保留线性项，忽略如果只保留线性项，忽略 h2 及以上各项，则：及以上各项，则：显式欧拉公式16第16页，本讲稿共65页局部截断误差的分析局部截断误差的分析利用泰勒公式展开，比较各算法与展开式的前几项利用泰勒公式展开，比较各算法与展开式的前几项将将 y(xn+1)在在 xn 点处用泰勒公式展开：点处用泰勒公式展开：显式欧拉法的局部截断误差：显式欧拉法的局部截断误差：欧拉法1 阶精度阶精度17第17页，本讲稿共65页补充：二元函数微分中值定理补充：二元函数微分中值定理18第18页
11、，本讲稿共65页y(xn+1)在 xn 点处展开：隐式欧拉法：1 阶精度阶精度19第19页，本讲稿共65页分别将分别将 y(xn+1),y(xn-1)在在 xn 点处用泰勒公式展开：点处用泰勒公式展开：二步欧拉法的局部截断误差：二步欧拉法的局部截断误差：二步欧拉法：2 阶精度阶精度20第20页，本讲稿共65页梯形公式欧拉法：y(xn+1)在 xn 点处展开：2 阶精度阶精度21第21页，本讲稿共65页各种欧拉法的比较各种欧拉法的比较方法精度评述显式欧拉法1最简单，精度低隐式欧拉法1不便计算，稳定性好二步欧拉法2需要两步初值，且第 2 个初值只能由其它方法给出，可能对后面的递推精度有影响梯形公式
12、欧 拉 法2精度有所提高，但为隐式，需要迭代求解，计算量大22第22页，本讲稿共65页改进的欧拉法改进的欧拉法从从上上述述例例子子可可以以看看到到，梯梯形形法法由由于于具具有有二二阶阶精精度度，其其局局部部截截断断误误差差比比显显式式欧欧拉拉法法和和隐隐式式欧欧拉拉法法小小，但但梯梯形形法实质上是一种隐式算法法实质上是一种隐式算法显显式式欧欧拉拉法法是是一一个个显显式式算算法法，虽虽然然计计算算量量较较小小，但但是是精度不高精度不高综综合合两两种种方方法法的的长长处处，可可以以先先用用显显式式欧欧拉拉法法求求出出 y(xn+1)的的一一个个粗粗略略近近似似值值，然然后后用用它它代代入入梯梯形形
13、法法公公式式的右端，用梯形法计算的右端，用梯形法计算 y(xn+1)的较为精确的近似值。的较为精确的近似值。23第23页，本讲稿共65页改进的欧拉法（续）改进的欧拉法（续）按照上述思想，可以建立如下预报按照上述思想，可以建立如下预报-校正系统：校正系统：按以上两式求解常微分方程的算法称为改进的欧拉按以上两式求解常微分方程的算法称为改进的欧拉法，它还可以表示为：法，它还可以表示为：嵌套形式平均化形式2 阶精度阶精度24第24页，本讲稿共65页用改进欧拉法求上例所述的初值问题并与欧拉法和用改进欧拉法求上例所述的初值问题并与欧拉法和梯形法比较误差的大小。梯形法比较误差的大小。解解：采用改进欧拉法的嵌
14、套形式：采用改进欧拉法的嵌套形式：25第25页，本讲稿共65页计算结果计算结果xn改进欧拉法 yn精确解 y(xn)误 差改进欧拉法欧拉法梯形法0.1 1.005000 1.004837 1.610-44.810-37.510-50.2 1.019205 1.019731 2.910-48.710-31.410-40.3 1.041218 1.040818 4.010-41.210-21.910-40.4 1.070802 1.070320 4.810-41.410-22.210-40.5 1.107076 1.106531 5.510-41.610-22.510-4可见，改进欧拉法的误差数量
15、级与梯形法大致相同，可见，改进欧拉法的误差数量级与梯形法大致相同，而比欧拉法小得多。而比欧拉法小得多。26第26页，本讲稿共65页改进的欧拉法的意义改进的欧拉法的意义改进的欧拉法的平均化形式y(xn+1)在点在点 xn 处的一阶展开式为：处的一阶展开式为：27第27页，本讲稿共65页改进的欧拉法的几何意义改进的欧拉法的几何意义P0 xyRhQ斜率:k2S斜率:k128第28页，本讲稿共65页龙格龙格-库塔库塔(Runge-Kutta)方法方法改进的欧拉法改进的欧拉法（2 2 阶精度）阶精度）y(xn+1)在点在点 xn 处的一阶泰勒展开式为：处的一阶泰勒展开式为：显式欧拉法显式欧拉法（1 1
16、阶精度）阶精度）29第29页，本讲稿共65页龙格龙格-库塔方法（续）库塔方法（续）显式欧拉法用一个点的值显式欧拉法用一个点的值 k1 作为作为 k*的近似值的近似值改改进进的的欧欧拉拉公公式式用用二二个个点点的的值值 k1 和和 k2 的的平平均均值值作作为为 k*近似值；近似值；改进的欧拉法比显式欧拉法精度高；改进的欧拉法比显式欧拉法精度高；在在 xn,xn+1 内内多多预预报报几几个个点点的的 ki 值值，并并用用其其加加权权平平均均值值作作为为 k*的的近近似似值值从从而而构构造造出出具具有有更更高高精精度度的的计计算算公公式，这就是式，这就是龙格龙格-库塔方法的基本思想。库塔方法的基本
17、思想。30第30页，本讲稿共65页二阶龙格二阶龙格-库塔方法库塔方法以以 k1 和和 k2 的加权平均来近似取代的加权平均来近似取代 k*为分析局部截断误差，令为分析局部截断误差，令 yn=y(xn)，由泰勒公式得：，由泰勒公式得：31第31页，本讲稿共65页补充：二元泰勒展开式补充：二元泰勒展开式32第32页，本讲稿共65页用二元泰勒公式展开将 k1,k2 代入 中可得：33第33页，本讲稿共65页二阶龙格二阶龙格-库塔方法（续）库塔方法（续）2 阶精度阶精度34第34页，本讲稿共65页四四个个未未知知变变量量，只只有有三三个方程，有无穷多组解个方程，有无穷多组解每每组组解解的的构构成成的的
18、龙龙格格-库库塔方法均为二阶塔方法均为二阶二二阶阶龙龙格格-库库塔塔方方法法即即为为改改 进进 的的 欧欧 拉拉 方方 法法变变形形的的欧欧拉拉法法中中 点点 方方 法法35第35页，本讲稿共65页三阶龙格三阶龙格-库塔方法库塔方法三三阶阶龙龙格格-库库塔塔方方法法是是用用三三个个值值 k1,k2,k3 的的加加权权平平均均来来近似取代近似取代 k*要要使使三三阶阶龙龙格格-库库塔塔方方法法具具有有三三阶阶精精度度，必必须须使使其其局局部部截断误差为截断误差为 O(h4)将将 k1,k2,k3 代代入入 yn+1 的的表表达达式式中中，在在(xn,yn)处处用用二二元元泰勒公式展开，与泰勒公式
19、展开，与 y(xn+1)在在 xn 处的泰勒展开式比较处的泰勒展开式比较36第36页，本讲稿共65页三阶龙格三阶龙格-库塔方法（续）库塔方法（续）类似二阶龙格类似二阶龙格-库塔方法的推导过程，库塔方法的推导过程，8 个待定系数个待定系数 c1,c2,c3,a2,a3,b21,b31,b32 应满足：应满足：8 个个未未知知参参数数，6 个个方方程程，有有无无穷穷多多组组解解库塔公式37第37页，本讲稿共65页四阶龙格四阶龙格-库塔方法库塔方法类似可以推出四阶龙格类似可以推出四阶龙格-库塔公式，常用的有：库塔公式，常用的有：标准四阶龙格-库塔公式38第38页，本讲稿共65页四阶龙格四阶龙格-库塔
20、方法（续）库塔方法（续）吉尔（Gill）公式4 阶阶以以上上龙龙格格-库库塔塔公公式式的的计计算算量量太太大大，并并且且精精度度不不一一定定提提高高，有有时时反反而而会会降降低低，因因此此实实际际应应用用中中一一般般选选用用四阶龙格四阶龙格-库塔已足可满足精度要求。库塔已足可满足精度要求。39第39页，本讲稿共65页用用经经典典四四阶阶龙龙格格-库库塔塔方方法法求求解解前前例例的的初初值值问问题题，并并与改进与改进 欧拉欧拉 法、梯形法在法、梯形法在 x5=0.5 处比较其误差大小处比较其误差大小解解：采用经典四阶龙格：采用经典四阶龙格-库塔库塔公式：公式：40第40页，本讲稿共65页四四阶阶
21、R-K方方法法的的精精度度比比二二阶阶方方法高得多法高得多精确解为：R-K方法的误差：改进欧拉法的误差：梯形法的误差：41第41页，本讲稿共65页变步长的龙格变步长的龙格-库塔方法库塔方法设 y(xn)在 xn 处的值 yn=y(xn)，当 xn+1=xn+h 时 y(xn+1)的近似值为 ，由于四阶 R-K 方法的精度为 4 阶，故局部截断误差为：用用四四阶阶R-K方方法法求求解解初初值值问问题题精精度度较较高高，但但要要从从理理论论上上给给出出误误差差|y(xn)-yn|的的估估计计式式则则比比较较困困难难；那那么么应应如如何判断计算结果的精度以及如何选择合适的步长何判断计算结果的精度以及
22、如何选择合适的步长 h？通通常常是是通通过过不不同同步步长长在在计计算算机机上上的的计计算算结结果果进进行行近近似似估计。估计。42第42页，本讲稿共65页若以 h/2 为步长，从 xn 出发，经过两步计算，得到y(xn+1)的近似值变步长的龙格变步长的龙格-库塔方法（续）库塔方法（续）以上每步的截断误差约为以上每步的截断误差约为 cn(h/2)5，于是两步的局部，于是两步的局部截断误差为：截断误差为：于是：整理得：43第43页，本讲稿共65页变步长的龙格变步长的龙格-库塔方法（续）库塔方法（续）记：，给定的精度要求为 eD e，反复将步长折半计算，直至 D e，取最终得到的 作为 y(xn+
23、1)的近似值。D D e e，再再将将步步长长折折半半一一次次计算，最终得到符合精度要求的计算，最终得到符合精度要求的 y(xn+1)的近似值。的近似值。44第44页，本讲稿共65页单步法的收敛性单步法的收敛性显式单步法可统一写成：显式单步法可统一写成：增增量量函函数数，仅仅依依赖赖于于函函数数 f，且仅仅是且仅仅是 xn,yn,h 的函数的函数求求 y =y(x)求 y(xn)，xn=x0+nh离散化求 yn y(xn)某种数值方法h 0 时，近似解时，近似解是否收敛到精确解是否收敛到精确解，它应当是一个固定节点，因此 h 0 时应同时附带 n 45第45页，本讲稿共65页单步法的收敛性（续
24、）单步法的收敛性（续）对于对于 p 阶的常微分方程数值算法，当阶的常微分方程数值算法，当 h 0，n 时，是否时，是否 yn+1 y(xn+1)？p 阶算法的局部截断误差为：显然：局部截断误差的前提假设是：局部截断误差局部截断误差 0 并不能保证算法收敛并不能保证算法收敛46第46页，本讲稿共65页单步法的收敛性（续）单步法的收敛性（续）定义：若求解某初值问题的单步数值法，对于固定的 当 h 0 且 n 时，它的近似 解趋向于精确解 y(xn)，即：则称该单步法是收敛的。则称该单步法是收敛的。定义定义：称：称 y(xn)-yn 为单步法的近似解为单步法的近似解 yn 的的整体截断整体截断 误差
25、误差。单步法收敛单步法收敛h 0 且 n 时，yn 的整体截断误差 0 47第47页，本讲稿共65页单步法的收敛性（续）单步法的收敛性（续）收敛性定理收敛性定理若某单步法满足以上条件，则该方法是收敛的若某单步法满足以上条件，则该方法是收敛的则该单步法的整体截断误差为：则该单步法的整体截断误差为：若单步法 具有 p 阶精度，且增量函数 关于 y 满足：Lipschitz 条件：初值初值 y0 是准确的是准确的48第48页，本讲稿共65页假设在前一步假设在前一步 yn 准确的前提下求得的近似值为：准确的前提下求得的近似值为：算法精度为算法精度为 p 阶，局部截断误差：阶，局部截断误差：49第49页
26、，本讲稿共65页50第50页，本讲稿共65页即：若初值是准确的，则若初值是准确的，则 e e 0=0，从而整体截断误差为：，从而整体截断误差为：y=e x 为单调增函数，当 时当 h 0 且 n 时，则51第51页，本讲稿共65页单步法的稳定性单步法的稳定性在讨论单步法收敛性时一般认为数值方法本身的计在讨论单步法收敛性时一般认为数值方法本身的计算过程是准确的，实际上并非如此：算过程是准确的，实际上并非如此：初始值初始值 y0 有误差有误差 =y0-y(x0)后续的每一步计算均有舍入误差后续的每一步计算均有舍入误差这些初始和舍入误差在计算过程的传播中是逐步衰这些初始和舍入误差在计算过程的传播中是
27、逐步衰减的还是恶性增长就是数值方法的稳定性问题减的还是恶性增长就是数值方法的稳定性问题52第52页，本讲稿共65页定义：若一种数值方法在节点 xn 处的数值解 yn 的扰动 ，而在以后各节点 ym(m n)上产生的扰动为 ，如果：单步法的稳定性（续）单步法的稳定性（续）定义：设在节点 xn 处用数值算法得到的理想数值解为 yn，而实际计算得到的近似解为 ，称差值：为第为第 n 步的数值解的步的数值解的扰动扰动。则称该数值方法是则称该数值方法是稳定稳定的。的。53第53页，本讲稿共65页单步法的稳定性（续）单步法的稳定性（续）欧拉法：由于函数 f (x,y)的多样性，数值稳定性的分析相当复杂，通
28、常只研究模型方程考察模型方程：即：假假设设在在节节点点值值 yn 上上有有扰扰动动 n，在在节节点点值值 yn+1 上上有有扰扰动动 n+1，且且 n+1 仅仅由由 n 引引起起（即即：计计算算过过程程中中不不再再引引起起新的误差）新的误差）54第54页，本讲稿共65页欧拉法稳定欧拉法稳定即：欧拉法稳定的条件：欧拉法稳定的条件：针对模型方程：的显式欧拉法：化简得：55第55页，本讲稿共65页隐式欧拉法：考察模型方程：即：化简为：假设 yn 上有扰动 ，则 yn+1 的扰动为：隐式欧拉法稳定隐式欧拉法稳定，上式均成立，所以：隐式欧拉法稳定是恒稳定的隐式欧拉法稳定是恒稳定的56第56页，本讲稿共6
29、5页一阶常微分方程组一阶常微分方程组57第57页，本讲稿共65页显式欧拉法显式欧拉法隐式欧拉法隐式欧拉法梯形公式欧拉法梯形公式欧拉法58第58页，本讲稿共65页改进的欧拉公式改进的欧拉公式59第59页，本讲稿共65页四阶龙格四阶龙格-库塔公式库塔公式60第60页，本讲稿共65页高阶微分方程的初值问题高阶微分方程的初值问题一般通过引入新的变量，将高阶微分方程化为一阶微一般通过引入新的变量，将高阶微分方程化为一阶微分方程组的方法进行求解分方程组的方法进行求解m 阶常微分方程：用变量替换可以将上述的高阶微分方程转化为一阶微分方程组用变量替换可以将上述的高阶微分方程转化为一阶微分方程组设：61第61页，本讲稿共65页初始条件为：则则 m 阶微分方程转化为如下的一阶微分方程转化为如下的一阶微分方程组：阶微分方程组：62第62页，本讲稿共65页转化例子转化例子将 转化为微分方程组。解：令 ，有：初始条件：初始条件：63第63页，本讲稿共65页64第64页，本讲稿共65页作业作业课本 P116：134691265第65页，本讲稿共65页