第4章常微分方程数值解法课件

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