常微分方程的数值解法ppt课件

第五章 常微分方程的数值解法,主要内容： 1、引言 2、欧拉方法 3、龙格库塔方法 4、单步法的收敛性和稳定性 5、线性多步法 6、一阶方程组与高阶方程,1,第一节 引言,在常微分方程课程里面讨论的是一些典型方程求解解析解的基本方法。 然而在生产实践和科学研究中遇到的微分方程往往比较复杂，在很多情况下，不能给出解的解析表达式；有时候即时能用解析表达式来表示，又因为计算量太大而不实用，有时候一些是已经有了求解的基本方法的典型方程，但实际使用时也是有困难的。 以上情况说明用求解解析解的基本方法来求微分方程的解往往是不适宜的，甚至很难办到。 实际问题中，对于求解微分方程，一般只要求得到解的若干个点上的近似值或者解的便于计算的近似表达式。 本章研究微分方程的数值解法，而且着重讨论微分方程中最简单的一类问题一阶方程的初值问题。,2,第一节 引言,1、一阶方程的初值问题 假定上式在区间a, b上存在唯一且足够 光滑的解y(x)。 所谓数值解法就是寻求解y(x)在一系列离散点，也称为节点处的值： 要计算出解函数 y(x) 在一系列节点 a = x0 x1 xn= b 处的近似值,3,第一节 引言,节点间距，即步长为： 通常采用等距节点，即hi = h (常数) 等间距节点 在这些节点上采用离散化方法（通常用数值积分、微分、泰勒展开等）将上述初值问题化成关于离散变量的相应问题。把这个相应问题的解yn作为y(xn)的近似值。这样求得的yn就是上述初值问题在节点xn上的数值解。一般说来，不同的离散化导致不同的方法。,4,第二节 欧拉方法,一、欧拉法Euler 1、向前差商近似导数,5,第二节 欧拉方法,2、举例 例1 用欧拉法求初值问题 当h = 0.02时在区间0, 0.10上的数值解。 解：根据欧拉公式可以得到： 此外，可以得到方程的真解：,6,第二节 欧拉方法,求解过程如下：,7,第二节 欧拉方法,3、欧拉方法的几何意义,根据已知条件：曲线y(x)上的点（x0,y0）及该点处曲线的导数f(x0,y0)，则可以得到过该点的直线：,该直线与xx1的交点P1，则P1的纵坐标y1为：,就用y1作为y(x1)的近似值,逐次进行后可以得到一条折线P0P1Pn，该折线看作是初值问题的积分曲线的近似，因此欧拉方法也称为欧拉折线法,8,第二节 欧拉方法,从上述几何意义上得知，由Euler法所得的折线明显偏离了积分曲线，可见此方法非常粗糙即误差太大。 4、欧拉法的局部截断误差 （1）截断误差定义 在假设 yi = y(xi)，即第 i 步计算是精确的前提下，考虑的截断误差 Ri+1 = y(xi+1) yi+1 ，称为局部截断误差 如图所示：APi+1即为欧拉方法在xi+1点的截断误差 （2）如果某种方法的局部截断误差是 则称该方法具有p阶精度,9,第二节 欧拉方法,（3）则截断误差的大小？ 写出y(xn+1)的泰勒展开式： 由欧拉方法可以得到： 则上面两个公式相减得到：,具有1阶精度,10,第二节 欧拉方法,二、改进的欧拉法 一阶方程的初值问题与如下积分方程是等价的： 当x = x1时 可以借助于数值积分，求y(x1)的值 1、用矩形公式,11,第二节 欧拉方法,可以推导出： 用矩形法计算右端的积分与用欧拉法计出的结果完全相同 2、用梯形公式 则可以推导出：,12,第二节 欧拉方法,梯形公式的截断误差：,梯形公式具有二阶精度，比欧拉方法有了进步,13,第二节 欧拉方法,和欧拉公式相比较，梯形公式在计算yi+1时候也只用到前一步的值yi，但是若yi已知，将yi带入公式求解时候，一般不能直接得到yi+1，而需要通过其他方法（比如迭代法）求解，所以梯形公式被称为隐式公式。 3、改进的欧拉方法 梯形公式是隐式的，一般用迭代法求解，计算量较大。实际中常将欧拉公式和梯形公式联合使用，先用欧拉公式得出一个y(xi+1)的近似值 称为预估值，然后对预估值使用梯形公式对它进行精确化，得到较为精确的近似值yi+1，称之为校正值，计算公式为： 这样的预估校正系统称为改进的欧拉方法。,14,第二节 欧拉方法,为了便于编写程序，常将上面的公式改写为如下式：,15,第二节 欧拉方法,4、举例 P90，例题5-1 在区间0, 1.5上，取h = 0.1，求解。 解：（1）用欧拉法计算公式如下： （2）用改进欧拉法计算公式如下：,16,本题的精确解为 ， 可用来检验近似解 的精确程度。 计算结果如表：,17,第二节 欧拉方法,P108 习题5-1,18,第二节 欧拉方法,5、欧拉两步公式 中心差商,19,第三节 龙格库塔方法,一、龙格库塔法的基本思想 1、平均斜率 考察差商： 根据微分中值定理：在闭区间a, b上连续，开区间(a, b)上可导，则至少存在(a, b)上的一点，使得下式成立： 根据上面公式可以得到：,称为区间xi,xi+1上的平均斜率,20,第三节 龙格库塔方法,因此只要对 K* 提供一种算法，就可以求得数值解，根据该观点对欧拉法及改进的欧拉进行分析。 2、基于平均斜率对欧拉法和改进的欧拉法进行分析 （1）在欧拉公式中，是取了一个点 xi上的斜率值f(xi,yi) 作为平均斜率K*的近似值的，已经知道其精度较低。 （2）对于改进的欧拉公式：,可以看出它用两个点xi 和 xi+1上的斜率K1和K2 的算术平均值作为平均斜率 K* 的近似值的，而xi+1处的斜率是由已知信息预测得到的。,21,第三节 龙格库塔方法,根据上面的分析得到：如果在区间 xi, xi+1上多测几个点的斜率值，然后取其加权平均作为平均斜率的近似值，有可能构造出具有更高精度的计算公式，此即为龙格库塔算法的基本思想。 二、二阶龙格库塔方法 1、推广改进的欧拉方法，考察区间xi, xi+1上的一点： 用xi 和xp 两个点上的斜率值K1和K2 的加权平均作为平均斜率 K*的近似值： 即取： 与改进的欧拉法类似，有： 如何得到xp的斜率？,22,第三节 龙格库塔方法,如何得到xp的斜率？ 根据改进的欧拉法，可以利用欧拉法预测 的值： 则可以得到点xi+p斜率K2： 则可以得到算法的具体表达式：,23,第三节 龙格库塔方法,2、选择参数使得算法具有2阶精度 计算上面公式的局部截断误差： 根据泰勒公式有： 另外：,24,25,第三节 龙格库塔方法,考虑到： 则有：,二元泰勒展开,两式相减可得到：,26,第三节 龙格库塔方法,两式相减可得到： 则要使得上式满足二阶精度，即Ri+1=O(h3)，只需要： 共有3个参数，但只需要满足2个条件，因此满足该式的参数不止一组，而是一簇，所有满足条件的公式通称为二阶龙格库塔公式。,27,第三节 龙格库塔方法,该公式的特别情况： 有： 此时，二阶龙格库塔公式就是改进的欧拉公式。,28,第三节 龙格库塔方法,如果p=1/2，则有： 此时二阶龙格库塔公式变成： 为变形的欧拉公式。,29,第三节 龙格库塔方法,3、高阶龙格库塔公式 为了进一步提高精度，可以考虑在区间xi, xi+1上，除了xi，xp外，再增加几个点。如再增加一个点：xi+m 利用该三点处的斜率的加权平均值作为K* 近似值： 此时计算公式为： 式中K1和K2的取法与前相同，如何获得K3？,30,第三节 龙格库塔方法,为了得到xi+m处的斜率K3，需要确定该点的函数值y(xi+m)，可以用二阶龙格库塔公式估计y(xi+m)的近似值yi+m： 然后通过计算函数值，就可以得到K3： 由此得到的 3阶龙格库塔 公式为：,31,第三节 龙格库塔方法,通过选择未知参数使得公式具有3阶精度，采用前面相同的方法，可以得到参数需满足如下条件：,32,第三节 龙格库塔方法,常用的三阶龙格库塔公式为：,其中,33,第三节 龙格库塔方法,经典的四阶龙格库塔公式为：,34,第三节 龙格库塔方法,35,