第四章常微分方程数值解法精选PPT.ppt
第四章常微分方程数值第四章常微分方程数值解法解法1 1第1页，此课件共65页哦常微分方程数值解法引言（常微分方程数值解法概述）引言（常微分方程数值解法概述）显式欧拉法、隐式欧拉法、二步欧拉法显式欧拉法、隐式欧拉法、二步欧拉法局部截断误差与精度局部截断误差与精度改进的欧拉方法改进的欧拉方法龙格龙格-库塔方法库塔方法收敛性与稳定性简述收敛性与稳定性简述一阶常微分方程组与高阶常微分方程一阶常微分方程组与高阶常微分方程2第2页，此课件共65页哦引言引言一阶常微分方程初值问题：一阶常微分方程初值问题：微分方程微分方程初始条件初始条件定理：定理：若若 f(x,y)在某闭区域在某闭区域 R：上上连连续续，且且在在 R 域域内内满满足足李李普普希希兹兹(Lipschitz)条条件件，即即存存在在正正数数 L，使使得得对对于于 R 域域内内的的任任意意两两值值 y1,y2，下下列列不不等等式式成立：成立：则上述初值问题的连续可微的解则上述初值问题的连续可微的解 y(x)存在并且唯一。存在并且唯一。3第3页，此课件共65页哦引言（续）引言（续）实实际际生生产产与与科科研研中中，除除少少数数简简单单情情况况能能获获得得初初值值问问题题的的初初等等解解（用用初初等等函函数数表表示示的的解解）外外，绝绝大大多多数数情情况况下是求不出初等解的。下是求不出初等解的。有有些些初初值值问问题题即即便便有有初初等等解解，也也往往往往由由于于形形式式过过于于复复杂而不便处理。杂而不便处理。实实用用的的方方法法是是在在计计算算机机上上进进行行数数值值求求解解：即即不不直直接接求求 y(x)的的显显式式解解，而而是是在在解解所所存存在在的的区区间间上上，求求得得一一系系列点列点 xn(n=0,1,2,)上解的近似值。上解的近似值。4第4页，此课件共65页哦欧拉欧拉（Euler）方法方法方法一方法一化导数为差商的方法化导数为差商的方法由于在逐步求解的过程中，由于在逐步求解的过程中，y(xn)的准确值无法求解出的准确值无法求解出来，因此用其近似值代替。来，因此用其近似值代替。为避免混淆，以下学习简记：为避免混淆，以下学习简记：y(xn)：待求函数：待求函数 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页哦欧拉方法（续）欧拉方法（续）方法二方法二数值积分法数值积分法同样以近似值同样以近似值 yn 代替精确代替精确值值 y(xn)可得可得：将微分方程将微分方程 y =f(x,y)在区在区间间 xn,xn+1 上上积积分：分：7第7页，此课件共65页哦欧拉方法的几何意义欧拉方法的几何意义xy08第8页，此课件共65页哦隐式欧拉法隐式欧拉法在数值积分法推导中，积分的近似值取为积分区间宽在数值积分法推导中，积分的近似值取为积分区间宽度与右端点处的函数值乘积，即：度与右端点处的函数值乘积，即：这样便得到了隐式欧拉法：这样便得到了隐式欧拉法：含有未知含有未知的函数值的函数值隐式欧拉法没有显式欧拉法方便隐式欧拉法没有显式欧拉法方便9第9页，此课件共65页哦二步欧拉法二步欧拉法在在数数值值积积分分法法推推导导中中，积积分分区区间间宽宽度度选选为为两两步步步步长长，即积分区间为：即积分区间为：xn-1,xn+1，则：，则：以以 y(x)在在 xn-1,xn 上的近似值上的近似值代替精确代替精确值值可得可得：需需要要前前两两步步的的计计算算结结果果中矩形公式中矩形公式10第10页，此课件共65页哦梯形公式欧拉法梯形公式欧拉法在数值积分法中，如果用梯形公式近似计算在数值积分法中，如果用梯形公式近似计算 f(x,y)在区间在区间 xn,xn+1 上的积分，即：上的积分，即：用近似值代替精确值可得梯形公式欧拉法：用近似值代替精确值可得梯形公式欧拉法：上式右端出现了未知项，可见梯形法是隐式欧拉法的一种；上式右端出现了未知项，可见梯形法是隐式欧拉法的一种；实际上，梯形公式欧拉法是显式欧拉法与隐式欧拉法的实际上，梯形公式欧拉法是显式欧拉法与隐式欧拉法的算术算术平均平均。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页哦例：（续）例：（续）隐式欧拉法：隐式欧拉法：化简得：化简得：梯形公式欧拉法：梯形公式欧拉法：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页哦局部截断误差局部截断误差为了简化分析某常微分方程数值算法的误差，现假为了简化分析某常微分方程数值算法的误差，现假设设 yn=y(xn)，即，即在前一步在前一步 yn 准确的前提下准确的前提下，估计：，估计：称上述误差称上述误差 Tn+1 为该常微分方程数值算法的为该常微分方程数值算法的局部截断误差局部截断误差如如果果某某个个常常微微分分方方程程数数值值算算法法的的局局部部截截断断误误差差可可表表示示为为 O(h p+1)，则称该数值算法的精度是，则称该数值算法的精度是 p 阶阶欧欧拉拉法法的的精精度度为为一一阶阶；二二步步欧欧拉拉法法的的精精度度为为二二阶阶；梯梯形形法的精度为二阶。法的精度为二阶。15第15页，此课件共65页哦泰勒展开法泰勒展开法如果初值问题中的如果初值问题中的 f(x,y)充分可微，则可将充分可微，则可将 y(xn+1)在点在点 xn 处展开：处展开：如果只保留线性项，忽略如果只保留线性项，忽略 h2 及以上各项，则：及以上各项，则：显式欧拉公式显式欧拉公式16第16页，此课件共65页哦局部截断误差的分析局部截断误差的分析利用泰勒公式展开，比较各算法与展开式的前几项利用泰勒公式展开，比较各算法与展开式的前几项将将 y(xn+1)在在 xn 点处用泰勒公式展开：点处用泰勒公式展开：显式欧拉法的局部截断误差：显式欧拉法的局部截断误差：欧拉法欧拉法1 阶精度阶精度17第17页，此课件共65页哦补充：二元函数微分中值定理补充：二元函数微分中值定理18第18页，此课件共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 个初值只能由其它方法给出，可能对后面的递推精度有影响梯形公式欧 拉 法2精度有所提高，但为隐式，需要迭代求解，计算量大22第22页，此课件共65页哦改进的欧拉法改进的欧拉法从从上上述述例例子子可可以以看看到到，梯梯形形法法由由于于具具有有二二阶阶精精度度，其其局局部部截截断断误误差差比比显显式式欧欧拉拉法法和和隐隐式式欧欧拉拉法法小小，但但梯梯形形法实质上是一种隐式算法法实质上是一种隐式算法显显式式欧欧拉拉法法是是一一个个显显式式算算法法，虽虽然然计计算算量量较较小小，但但是是精度不高精度不高综综合合两两种种方方法法的的长长处处，可可以以先先用用显显式式欧欧拉拉法法求求出出 y(xn+1)的的一一个个粗粗略略近近似似值值，然然后后用用它它代代入入梯梯形形法法公公式式的右端，用梯形法计算的右端，用梯形法计算 y(xn+1)的较为精确的近似值。的较为精确的近似值。23第23页，此课件共65页哦改进的欧拉法（续）改进的欧拉法（续）按照上述思想，可以建立如下预报按照上述思想，可以建立如下预报-校正系统：校正系统：按以上两式求解常微分方程的算法称为改进的欧拉法，它按以上两式求解常微分方程的算法称为改进的欧拉法，它还可以表示为：还可以表示为：嵌套形式嵌套形式平均化形式平均化形式2 阶精度阶精度24第24页，此课件共65页哦用改进欧拉法求上例所述的初值问题并与欧拉法和用改进欧拉法求上例所述的初值问题并与欧拉法和梯形法比较误差的大小。梯形法比较误差的大小。解解：采用改进欧拉法的嵌套形式：采用改进欧拉法的嵌套形式：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可见，改进欧拉法的误差数量级与梯形法大致相同，而比欧拉可见，改进欧拉法的误差数量级与梯形法大致相同，而比欧拉法小得多。法小得多。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 阶精度）阶精度）29第29页，此课件共65页哦龙格龙格-库塔方法（续）库塔方法（续）显式欧拉法用一个点的值显式欧拉法用一个点的值 k1 作为作为 k*的近似值的近似值改改进进的的欧欧拉拉公公式式用用二二个个点点的的值值 k1 和和 k2 的的平平均均值值作作为为 k*近似值；近似值；改进的欧拉法比显式欧拉法精度高；改进的欧拉法比显式欧拉法精度高；在在 xn,xn+1 内内多多预预报报几几个个点点的的 ki 值值，并并用用其其加加权权平平均均值值作作为为 k*的的近近似似值值从从而而构构造造出出具具有有更更高高精精度度的的计计算算公公式，这就是式，这就是龙格龙格-库塔方法的基本思想。库塔方法的基本思想。30第30页，此课件共65页哦二阶龙格二阶龙格-库塔方法库塔方法以以 k1 和和 k2 的加权平均来近似取代的加权平均来近似取代 k*为分析局部截断误差，令为分析局部截断误差，令 yn=y(xn)，由泰勒公式得：，由泰勒公式得：31第31页，此课件共65页哦补充：二元泰勒展开式补充：二元泰勒展开式32第32页，此课件共65页哦用二元泰勒公式展开用二元泰勒公式展开将将 k1,k2 代入代入 中可得：中可得：33第33页，此课件共65页哦二阶龙格二阶龙格-库塔方法（续）库塔方法（续）2 阶精度阶精度34第34页，此课件共65页哦四四个个未未知知变变量量，只只有有三三个方程，有无穷多组解个方程，有无穷多组解每每组组解解的的构构成成的的龙龙格格-库库塔方法均为二阶塔方法均为二阶二二阶阶龙龙格格-库库塔塔方方法法即即为为 改改 进进 的的 欧欧 拉拉 方方 法法变变形形的的欧欧拉拉法法中中 点点 方方 法法35第35页，此课件共65页哦三阶龙格三阶龙格-库塔方法库塔方法三三阶阶龙龙格格-库库塔塔方方法法是是用用三三个个值值 k1,k2,k3 的的加加权权平平均均来来近似取代近似取代 k*要要使使三三阶阶龙龙格格-库库塔塔方方法法具具有有三三阶阶精精度度，必必须须使使其其局局部部截截断误差为断误差为 O(h4)将将 k1,k2,k3 代代入入 yn+1 的的表表达达式式中中，在在(xn,yn)处处用用二二元元泰泰勒勒公公式展开，与式展开，与 y(xn+1)在在 xn 处的泰勒展开式比较处的泰勒展开式比较36第36页，此课件共65页哦三阶龙格三阶龙格-库塔方法（续）库塔方法（续）类似二阶龙格类似二阶龙格-库塔方法的推导过程，库塔方法的推导过程，8 个待定系数个待定系数 c1,c2,c3,a2,a3,b21,b31,b32 应满足：应满足：8 个个未未知知参参数数，6 个个方方程程，有有 无无 穷穷 多多 组组 解解库塔公式库塔公式37第37页，此课件共65页哦四阶龙格四阶龙格-库塔方法库塔方法类似可以推出四阶龙格类似可以推出四阶龙格-库塔公式，常用的有：库塔公式，常用的有：标准四阶龙格标准四阶龙格-库塔公式库塔公式38第38页，此课件共65页哦四阶龙格四阶龙格-库塔方法（续）库塔方法（续）吉尔（吉尔（Gill）公式）公式4 阶阶以以上上龙龙格格-库库塔塔公公式式的的计计算算量量太太大大，并并且且精精度度不不一一定定提提高高，有有时时反反而而会会降降低低，因因此此实实际际应应用用中中一一般般选选用用四四阶阶龙龙格格-库塔已足可满足精度要求。库塔已足可满足精度要求。39第39页，此课件共65页哦用用经经典典四四阶阶龙龙格格-库库塔塔方方法法求求解解前前例例的的初初值值问问题题，并并与改进与改进 欧拉欧拉 法、梯形法在法、梯形法在 x5=0.5 处比较其误差大小处比较其误差大小解解：采用经典四阶龙格：采用经典四阶龙格-库塔库塔公式：公式：40第40页，此课件共65页哦四四阶阶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|的的估估计计式式则则比比较较困困难难；那那么么应应如如何何判判断计算结果的精度以及如何选择合适的步长断计算结果的精度以及如何选择合适的步长 h？通常是通过不同步长在计算机上的计算结果进行近似估计。通常是通过不同步长在计算机上的计算结果进行近似估计。42第42页，此课件共65页哦若以若以 h/2 为步长，从为步长，从 xn 出发，经过两步计算，得到出发，经过两步计算，得到y(xn+1)的近似值的近似值变步长的龙格变步长的龙格-库塔方法（续）库塔方法（续）以上每步的截断误差约为以上每步的截断误差约为 cn(h/2)5，于是两步的局部截断，于是两步的局部截断误差为：误差为：于是：于是：整理得：整理得：43第43页，此课件共65页哦变步长的龙格变步长的龙格-库塔方法（续）库塔方法（续）记：记：，给定的精度要求为，给定的精度要求为 e eD D e e，反反复复将将步步长长折折半半计计算算，直直至至 D D e e，取取最最终终得得到的到的 作为作为 y(xn+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页哦单步法的收敛性（续）单步法的收敛性（续）对于对于 p 阶的常微分方程数值算法，当阶的常微分方程数值算法，当 h 0，n 时，是否时，是否 yn+1 y(xn+1)？p 阶算法的局部截断误差为：阶算法的局部截断误差为：显然：显然：局部截断误差的前提假设是：局部截断误差的前提假设是：局部截断误差局部截断误差 0 并不能保证算法收敛并不能保证算法收敛46第46页，此课件共65页哦单步法的收敛性（续）单步法的收敛性（续）定义定义：若求解某初值问题的单步数值法，对于固定的：若求解某初值问题的单步数值法，对于固定的 当当 h 0 且且 n 时，它的近似时，它的近似 解趋向于精确解解趋向于精确解 y(xn)，即：，即：则称该单步法是收敛的。则称该单步法是收敛的。定义定义：称：称 y(xn)-yn 为单步法的近似解为单步法的近似解 yn 的的整体截断整体截断 误差误差。单步法收敛单步法收敛h 0 且且 n 时，时，yn 的整体截断误差的整体截断误差 0 47第47页，此课件共65页哦单步法的收敛性（续）单步法的收敛性（续）收敛性定理收敛性定理若某单步法满足以上条件，则该方法是收敛的若某单步法满足以上条件，则该方法是收敛的则该单步法的整体截断误差为：则该单步法的整体截断误差为：若若单单步步法法 具具有有 p 阶阶精精度度，且增量函数且增量函数 关于关于 y 满足：满足：Lipschitz 条件：条件：初值初值 y0 是准确的是准确的48第48页，此课件共65页哦假设在前一步假设在前一步 yn 准确的前提下求得的近似值为：准确的前提下求得的近似值为：算法精度为算法精度为 p 阶，局部截断误差：阶，局部截断误差：49第49页，此课件共65页哦50第50页，此课件共65页哦即：即：若初值是准确的，则若初值是准确的，则 e e 0=0，从而整体截断误差为：，从而整体截断误差为：y=e x 为单调增函数，当为单调增函数，当 时时当当 h 0 且且 n 时时，则，则51第51页，此课件共65页哦单步法的稳定性单步法的稳定性在讨论单步法收敛性时一般认为数值方法本身的计在讨论单步法收敛性时一般认为数值方法本身的计算过程是准确的，实际上并非如此：算过程是准确的，实际上并非如此：初始值初始值 y0 有误差有误差 =y0-y(x0)后续的每一步计算均有舍入误差后续的每一步计算均有舍入误差这些初始和舍入误差在计算过程的传播中是逐步衰减这些初始和舍入误差在计算过程的传播中是逐步衰减的还是恶性增长就是数值方法的稳定性问题的还是恶性增长就是数值方法的稳定性问题52第52页，此课件共65页哦定定义义：若若一一种种数数值值方方法法在在节节点点 xn 处处的的数数值值解解 yn 的的扰扰动动 ，而而在在以以后后各各节节点点 ym(m n)上上产产生生的的扰扰动动为为 ，如果：，如果：单步法的稳定性（续）单步法的稳定性（续）定义定义：设在节点：设在节点 xn 处用数值算法得到的理想数值解为处用数值算法得到的理想数值解为 yn，而实际计算得到的近似解为，而实际计算得到的近似解为 ，称差值：，称差值：为第为第 n 步的数值解的步的数值解的扰动扰动。则称该数值方法是则称该数值方法是稳定稳定的。的。53第53页，此课件共65页哦单步法的稳定性（续）单步法的稳定性（续）欧拉法：欧拉法：由于函数由于函数 f (x,y)的多样性，数值稳定性的分析相当复的多样性，数值稳定性的分析相当复杂，通常只研究模型方程杂，通常只研究模型方程考察模型方程：考察模型方程：即：即：假假设设在在节节点点值值 yn 上上有有扰扰动动 n，在在节节点点值值 yn+1 上上有有扰扰动动 n+1，且且 n+1 仅仅由由 n 引引起起（即即：计计算算过过程程中中不不再再引引起起新新的误差）的误差）54第54页，此课件共65页哦欧拉法稳定欧拉法稳定即：即：欧拉法稳定的条件：欧拉法稳定的条件：针对模型方程：针对模型方程：的显式欧拉法：的显式欧拉法：化简得：化简得：55第55页，此课件共65页哦隐式欧拉法：隐式欧拉法：考察模型方程：考察模型方程：即：即：化简为：化简为：假设假设 yn 上有扰动上有扰动 ，则，则 yn+1 的扰动为：的扰动为：隐式欧拉法稳定隐式欧拉法稳定，上式均成立，所以：，上式均成立，所以：隐式欧拉法稳定是恒稳定的隐式欧拉法稳定是恒稳定的56第56页，此课件共65页哦一阶常微分方程组一阶常微分方程组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页哦
