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

1、数值分析Numerical Analysis第八章常微分方程数值解法郑州大学研究生课程 （2010-2011学年第一学期） ISCM 2007,Beijing China1第八章 常微分方程数值解法 8.1 引言8.2 欧拉(Euler)法8.3 改进欧拉(Euler)方法8.4 单步法的稳定性ISCM 2007,Beijing China/69 郑州大学研究生2010-2011学年课程 数值分析 Numerical Analysis8.1 引言问题提出 倒葫芦形状容器壁上的刻度问题.对于圆柱形状容器壁上的容积刻度,可以利用圆柱体体积公式其中直径D为常数.由于体积V与相对于容器底部的任意高度H

2、的函数关系明确,因此在容器上可以方便地标出容器刻度,而对于几何形状不是规则的容器,比如倒葫芦形状容器壁上如何标出刻度呢?ISCM 2007,Beijing China/69 郑州大学研究生2010-2011学年课程 数值分析 Numerical Analysis8.1 引言下表是经过测量得到部分容器高度与直径的关系.H 0 0.2 0.4 0.6 0.8 1.0D 0 0.11 0.26 0.56 1.04 1.17根据上表的数据,可以拟合出倒葫芦形状容器的图,建立如图所示的坐标轴后,问题即为如何根据任意高度x标出容器体积V的刻度,由微元思想分析可知ISCM 2007,Beijing Chin

3、a/69 郑州大学研究生2010-2011学年课程 数值分析 Numerical Analysis8.1 引言其中x表示高度,直径D是高度x的函数,记为D(x),因此得到如下微分方程初值问题只要求解上述方程,就可求出体积V与高度x之间的函数关系,从而可标出容器壁上容积的刻度,但问题是函数D(x)无解析表达式,我们无法求出其解析解.ISCM 2007,Beijing China/69 郑州大学研究生2010-2011学年课程 数值分析 Numerical Analysis8.1 引言 包含自变量、未知函数及未知函数的导数或微分的方程称为微分方程。在微分方程中, 自变量的个数只有一个, 称为常微分

4、方程。自变量的个数为两个或两个以上的微分方程叫偏微分方程。微分方程中出现的未知函数最高阶导数的阶数称为微分方程的阶数。如果未知函数y及其各阶导数都是一次的,则称它是线性的,否则称为非线性的。 ISCM 2007,Beijing China/69 郑州大学研究生2010-2011学年课程 数值分析 Numerical Analysis常微分方程( ODEs 未知函数是一元函数) 偏微分方程( PDEs 未知函数是多元函数) ISCM 2007,Beijing China/69 郑州大学研究生2010-2011学年课程 数值分析 Numerical Analysis同一个微分方程,具有不同的初始条

5、件ISCM 2007,Beijing China/69 郑州大学研究生2010-2011学年课程 数值分析 Numerical Analysis当x=0时,y=1,可得c=1特解当x=0时,y=1,可得c=-1特解两边积分通解ISCM 2007,Beijing China/69 郑州大学研究生2010-2011学年课程 数值分析 Numerical Analysis8.1 引言 在高等数学中，对于常微分方程的求解，给出了一些典型方程求解析解的基本方法，如可分离变量法、常系数齐次线性方程的解法、常系数非齐次线性方程的解法等。但能求解的常微分方程仍然是有限的，大多数的常微分方程是不可能给出解析解。

6、 ISCM 2007,Beijing China/69 郑州大学研究生2010-2011学年课程 数值分析 Numerical Analysis8.1 引言 待求解的问题：一阶常微分方程的初值问题 /* Initial-Value Problem */: 解的存在唯一性（“常微分方程”理论）：只要 f (x, y) 在a, b R1 上连续，且关于 y 满足 Lipschitz 条件，即存在与 x, y 无关的常数 L 使对任意定义在 a, b 上的 y1(x) 和 y2(x) 都成立，则上述IVP存在唯一解。ISCM 2007,Beijing China/69 郑州大学研究生2010-201

7、1学年课程 数值分析 Numerical Analysis解析解法：（常微分方程理论）只能求解极少一类常微分方程；实际中给定的问题不一定是解析表达式，而是函数表，无法用解析解法。如何求解ISCM 2007,Beijing China/69 郑州大学研究生2010-2011学年课程 数值分析 Numerical AnalysisISCM 2007,Beijing China/69 郑州大学研究生2010-2011学年课程 数值分析 Numerical AnalysisISCM 2007,Beijing China/69 郑州大学研究生2010-2011学年课程 数值分析 Numerical An

8、alysis8.2 欧拉(Euler)法 欧拉（Euler）方法是解初值问题的最简单的数值方法。初值问题的解y=y(x)代表通过点 的一条称之为微分方程的积分曲线。积分曲线上每一点 的切线的斜率 等于函数 在这点的值。 ISCM 2007,Beijing China/69 郑州大学研究生2010-2011学年课程 数值分析 Numerical AnalysisEuler法的求解过程是:从初始点P0(即点(x0,y0)出发,作积分曲线y=y(x)在P0点上切线 (其斜率为 ),与x=x1直线相交于P1点(即点(x1,y1),得到y1作为y(x1)的近似值,如上图所示。过点(x0,y0),以f(x

9、0,y0)为斜率的切线方程为 当 时,得 这样就获得了P1点的坐标。 ISCM 2007,Beijing China/69 郑州大学研究生2010-2011学年课程 数值分析 Numerical Analysis同样, 过点P1(x1,y1),作积分曲线y=y(x)的切线交直线x=x2于P2点,切线 的斜率 直线方程为当 时,得 ISCM 2007,Beijing China/69 郑州大学研究生2010-2011学年课程 数值分析 Numerical Analysis当 时,得由此获得了P2的坐标。重复以上过程,就可获得一系列的点:P1,P1,Pn。对已求得点以 为斜率作直线 取ISCM 2

10、007,Beijing China/69 郑州大学研究生2010-2011学年课程 数值分析 Numerical Analysis 从图形上看,就获得了一条近似于曲线y=y(x) 的折线 。这样,从x0逐个算出对应的数值解 ISCM 2007,Beijing China/69 郑州大学研究生2010-2011学年课程 数值分析 Numerical Analysis8.2 欧拉(Euler)法通常取 (常数),则Euler法的计算格式 i=0,1,n ( 8.2 ) ISCM 2007,Beijing China/69 郑州大学研究生2010-2011学年课程 数值分析 Numerical An

11、alysis8.2 欧拉(Euler)法还可用以下方法推导Euler格式： 数值微分 数值积分法对微分方程的离散，可以有多种思路，但最基本的想法是“以直代曲”ISCM 2007,Beijing China/69 郑州大学研究生2010-2011学年课程 数值分析 Numerical Analysis8.2 欧拉(Euler)法(1) 用差商近似导数差分方程初值问题向前Euler方法ISCM 2007,Beijing China/69 郑州大学研究生2010-2011学年课程 数值分析 Numerical Analysis8.2 欧拉(Euler)法若用向后差商近似导数，即向后Euler方法IS

12、CM 2007,Beijing China/69 郑州大学研究生2010-2011学年课程 数值分析 Numerical Analysis8.2 欧拉(Euler)法（2）用数值积分方法ISCM 2007,Beijing China/69 郑州大学研究生2010-2011学年课程 数值分析 Numerical Analysis8.2 欧拉(Euler)法若对积分用梯形公式，则得梯形欧拉公式ISCM 2007,Beijing China/69 郑州大学研究生2010-2011学年课程 数值分析 Numerical Analysis例8.2.1 用欧拉法解初值问题 取步长h=0.2 ,计算过程保留

13、4位小数 解: h=0.2, 欧拉迭代格式 当 k=0, x1=0.2时，已知x0=0，y0=1，有 y(0.2)y1=0.21(401)0.8当 k=1, x2=0.4时，已知x1 =0.2， y1 =0.8，有 y(0.4) y2 =0.20.8(40.20.8)0.6144当 k=2, x3 =0.6时，已知x2 =0.4， y2 =0.6144，有 y(0.6) y3=0.20.6144(4-0.40.6144)=0.4613 ISCM 2007,Beijing China/69 郑州大学研究生2010-2011学年课程 数值分析 Numerical Analysis8.2 欧拉(Eu

14、ler)法的解作为微分方程初值问题的数值解，即以差分方程初值问题ISCM 2007,Beijing China/69 郑州大学研究生2010-2011学年课程 数值分析 Numerical Analysis8.2 欧拉(Euler)法ISCM 2007,Beijing China/69 郑州大学研究生2010-2011学年课程 数值分析 Numerical Analysis8.2 欧拉(Euler)法x0 x1x2x3y0h h h 欧拉折线法ISCM 2007,Beijing China/69 郑州大学研究生2010-2011学年课程 数值分析 Numerical Analysis解：Eul

15、er公式为当h=0.5时ISCM 2007,Beijing China/69 郑州大学研究生2010-2011学年课程 数值分析 Numerical Analysis当h=0.25时ISCM 2007,Beijing China/69 郑州大学研究生2010-2011学年课程 数值分析 Numerical Analysis00.50.751.010.25h = 0.5h = 0.25ISCM 2007,Beijing China/69 郑州大学研究生2010-2011学年课程 数值分析 Numerical Analysis8.2 欧拉(Euler)法欧拉方法的收敛性ISCM 2007,Beij

16、ing China/69 郑州大学研究生2010-2011学年课程 数值分析 Numerical Analysis8.2 欧拉(Euler)法局部截断误差称为局部截断误差ISCM 2007,Beijing China/69 郑州大学研究生2010-2011学年课程 数值分析 Numerical Analysis8.2 欧拉(Euler)法欧拉方法的收敛性定义 若给定方法的局部截断误差满足则称该方法是 P 阶的，或称为具有 P 阶精度。ISCM 2007,Beijing China/69 郑州大学研究生2010-2011学年课程 数值分析 Numerical Analysis8.2 欧拉(Eul

17、er)法整体截断误差ISCM 2007,Beijing China/69 郑州大学研究生2010-2011学年课程 数值分析 Numerical Analysis8.2 欧拉(Euler)法欧拉方法的收敛性ISCM 2007,Beijing China/69 郑州大学研究生2010-2011学年课程 数值分析 Numerical Analysis由此知，当 ISCM 2007,Beijing China/69 郑州大学研究生2010-2011学年课程 数值分析 Numerical Analysis8.2 欧拉(Euler)法 注 整体截断误差与局部截断误差的关系： ISCM 2007,Beij

18、ing China/69 郑州大学研究生2010-2011学年课程 数值分析 Numerical Analysis8.2 欧拉(Euler)法 向后欧拉公式隐式欧拉法或向后欧拉法 /* implicit Euler method or backward Euler method*/xn+1点向后差商近似导数隐式或后退欧拉公式ISCM 2007,Beijing China/69 郑州大学研究生2010-2011学年课程 数值分析 Numerical Analysis8.2 欧拉(Euler)法 向后欧拉公式由于未知数 yn+1 同时出现在等式的两边，故称为隐式 /* implicit */ 欧拉

19、公式，而前者称为显式 /* explicit */ 欧拉公式。隐式公式不能直接求解，一般需要用Euler显式公式得到初值，然后用Euler隐式公式迭代求解。因此隐式公式较显式公式计算复杂，但稳定性好（后面分析）。 隐式欧拉公式中的未知数 yn+1 可通过以下迭代法求解：ISCM 2007,Beijing China/69 郑州大学研究生2010-2011学年课程 数值分析 Numerical Analysis8.2 欧拉(Euler)法 向后欧拉公式迭代法求隐式欧拉格式中yn+1的收敛性ISCM 2007,Beijing China/69 郑州大学研究生2010-2011学年课程 数值分析 N

20、umerical Analysis 见上图， 显然，这种近似也有一定误差，如何估计这种误差y(xn+1) yn+1 ？方法同上，基于Taylor展开估计局部截断误差。但是注意，隐式公式中右边含有f(xn+1 , yn +1 ) ,由于yn +1不准确，所以不能直接用y (xn+1)代替f(xn+1 , yn +1 ) 设已知曲线上一点 Pn (xn , yn ),过该点作弦线，斜率为(xn+1 , yn +1 ) 点的方向场f(x,y)。若步长h充分小，可用弦线和垂线x=xn+1的交点近似曲线与垂线的交点。几何意义xnxn+1PnPn+1xyy(x)ISCM 2007,Beijing Chin

21、a/69 郑州大学研究生2010-2011学年课程 数值分析 Numerical Analysis 隐式欧拉法的局部截断误差：ISCM 2007,Beijing China/69 郑州大学研究生2010-2011学年课程 数值分析 Numerical AnalysisISCM 2007,Beijing China/69 郑州大学研究生2010-2011学年课程 数值分析 Numerical Analysis 隐式欧拉法的局部截断误差：即隐式欧拉公式具有 1 阶精度。 隐式欧拉法的局部截断误差：ISCM 2007,Beijing China/69 郑州大学研究生2010-2011学年课程 数值分

22、析 Numerical Analysis8.2 欧拉(Euler)法 向后欧拉公式比较欧拉显式公式和隐式公式及其局部截断误差显式公式隐式公式ISCM 2007,Beijing China/69 郑州大学研究生2010-2011学年课程 数值分析 Numerical Analysis 若将这两种方法进行算术平均，即可消除误差的主要部分/*leading term*/而获得更高的精度,称为梯形法 梯形公式 /* trapezoid formula */ 显、隐式两种算法的平均注：的确有局部截断误差 ， 即梯形公式具有2 阶精度，比欧拉方法有了进步。但注意到该公式是隐式公式，计算时不得不用到迭代法，

23、其迭代收敛性与欧拉公式相似。ISCM 2007,Beijing China/69 郑州大学研究生2010-2011学年课程 数值分析 Numerical Analysis例8.2.3 对初值问题 证明用梯形公式求得的近似解为 并证明当步长h0时,yn收敛于精确解证明: 解初值问题的梯形公式为 整理成显式 反复迭代,得到 ISCM 2007,Beijing China/69 郑州大学研究生2010-2011学年课程 数值分析 Numerical Analysis公式局部截断误差精度显隐稳定性步数欧拉显式公式1阶显差单步欧拉隐式公式1阶隐好单步梯形公式2阶隐好单步欧拉法小结ISCM 2007,Be

24、ijing China/69 郑州大学研究生2010-2011学年课程 数值分析 Numerical Analysis8.3 改进欧拉(Euler)方法ISCM 2007,Beijing China/69 郑州大学研究生2010-2011学年课程 数值分析 Numerical Analysis8.3 改进欧拉(Euler)方法ISCM 2007,Beijing China/69 郑州大学研究生2010-2011学年课程 数值分析 Numerical Analysis8.3 改进欧拉(Euler)方法 显式欧拉公式计算工作量小,但精度低。梯形公式虽提高了精度,但为隐式公式,需用迭代法求解,计算工

25、作量大。综合欧拉公式和梯形公式便可得到改进的欧拉公式。 结合已有格式的优点，以得到计算方便、计算量减少且精度保持的数值格式ISCM 2007,Beijing China/69 郑州大学研究生2010-2011学年课程 数值分析 Numerical Analysis8.3 改进欧拉(Euler)方法 先用欧拉公式(8.2)求出一个初步的近似值,称为预测值, 它的精度不高, 再用梯形公式对它校正一次,即迭代一次,求得yn+1,称为校正值, 这种预测-校正方法称为改进的欧拉公式：称为Euler公式与梯形公式的预测校正系统。ISCM 2007,Beijing China/69 郑州大学研究生2010-

26、2011学年课程 数值分析 Numerical Analysis8.3 改进欧拉(Euler)方法实际计算时，常改写成以下形式几何解释xnxn+1ABPn+1=(A+B)/2欧拉法改进欧拉法梯形法ISCM 2007,Beijing China/69 郑州大学研究生2010-2011学年课程 数值分析 Numerical AnalysispredictorcorrectorISCM 2007,Beijing China/69 郑州大学研究生2010-2011学年课程 数值分析 Numerical Analysis8.3 改进欧拉(Euler)方法 可以证明,改进的欧拉公式的精度为二阶。这是一种一

27、步显式格式,它可以表示为嵌套形式。ISCM 2007,Beijing China/69 郑州大学研究生2010-2011学年课程 数值分析 Numerical Analysis例8.3.1ISCM 2007,Beijing China/69 郑州大学研究生2010-2011学年课程 数值分析 Numerical AnalysisISCM 2007,Beijing China/69 郑州大学研究生2010-2011学年课程 数值分析 Numerical Analysis8.3 改进欧拉(Euler)方法ISCM 2007,Beijing China/69 郑州大学研究生2010-2011学年课程

28、 数值分析 Numerical Analysis改进欧拉法的算法ISCM 2007,Beijing China/69 郑州大学研究生2010-2011学年课程 数值分析 Numerical Analysis8.4 单步法的稳定性 稳定性在微分方程的数值解法中是一个非常重要的问题。因为微分方程初值问题的数值方法是用差分格式进行计算的，而在差分方程的求解过程中，存在着各种计算误差，这些计算误差如舍入误差等引起的扰动，在传播过程中，可能会大量积累，对计算结果的准确性将产生影响。这就涉及到算法稳定性问题。 ISCM 2007,Beijing China/69 郑州大学研究生2010-2011学年课程 数值分析 Numerical Analysis例：考察初值问题 在区间0, 0.5上的解。分别用欧拉显、隐式格式和改进的欧拉格式计算数值解。0.00.10.20.30.40.5精确解改进欧拉法 欧拉隐式欧拉显式 节点 xi 1.00002.0000 4.00008.0000 1.6000101 3.2000101 1.00002.5000101 6.25001021.56251023.90631039.76561041.00002.50006.25001.56261013.9