常微分方程初值问题数值解法精.ppt

《常微分方程初值问题数值解法精.ppt》由会员分享，可在线阅读，更多相关《常微分方程初值问题数值解法精.ppt（69页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、常微分方程初值问题数值解法第1页，本讲稿共69页第九章第九章 常微分方程初值问题常微分方程初值问题数值解法数值解法/*Numerical Method for Ordinary Differential Equations*/许多实际问题的数学模型是微分方程或微分方程许多实际问题的数学模型是微分方程或微分方程的初值问题的初值问题,如物体运动如物体运动,电路震荡电路震荡,化学反映及生物群化学反映及生物群体的变化等。体的变化等。能用解析方法求出精确解的微分方程为数不多能用解析方法求出精确解的微分方程为数不多,而且而且有的方程即使有解析解有的方程即使有解析解,也可能由于解的表达式非常复杂也可能由于解

2、的表达式非常复杂而不易计算而不易计算,因此有必要研究微分方程的数值解法。因此有必要研究微分方程的数值解法。第2页，本讲稿共69页常微分方程中介绍的微分方程主要有常微分方程中介绍的微分方程主要有:(1)(1)变量可分离的方程变量可分离的方程 (2)(2)一阶线性微分方程一阶线性微分方程(贝努利方程贝努利方程)(3)(3)可降阶的一类高阶方程可降阶的一类高阶方程 (4)(4)二阶常系数齐次微分方程二阶常系数齐次微分方程 (5)(5)二阶常系数非齐次微分方程二阶常系数非齐次微分方程 (6)(6)全微分方程全微分方程本章主要介绍本章主要介绍一阶一阶常微分方程初值问题的常微分方程初值问题的数值解法数值解

3、法。第3页，本讲稿共69页图形解 xyo简单的微分方程复杂、大型的微分方程解析解 y=f(x)数值解(xi,yi)欧拉方法改进欧拉方法 梯形法龙格-库塔法第4页，本讲稿共69页 初值问题及其初值问题及其数值解数值解的概念的概念1 引言引言常用的一些常用的一些解析解法解析解法：常数常数变易变易法、法、Lapalace变换等变换等分离变量分离变量法、变量法、变量代换代换、一阶一阶常微分方程初值问题：常微分方程初值问题：第5页，本讲稿共69页对于初值问题对于初值问题 ，如果，如果 在下列区域内连续：在下列区域内连续：（解的（解的存在唯一存在唯一性）性）且关于且关于 满足满足Lipschitz条件，即

4、存在常数条件，即存在常数 ，使，使则初值问题则初值问题 存在唯一解，且解是存在唯一解，且解是连续可微连续可微的。的。所谓所谓数值解数值解是指：在解的是指：在解的存在区间存在区间上取一系列点上取一系列点逐个求出逐个求出 的近似值的近似值等距等距节点：节点：步长步长第6页，本讲稿共69页 初值问题初值问题 的解析解及其数值解的的解析解及其数值解的几何几何意义：意义：初值问题初值问题 的解表示过点的解表示过点 的一条的一条曲线曲线初值问题初值问题 的数值解表示一组的数值解表示一组离散点列离散点列可用可用拟合拟合方法求该组数据方法求该组数据 的的近似曲线近似曲线积分积分曲线曲线第7页，本讲稿共69页建

5、立微分方程数值解法建立微分方程数值解法,首先要将微分方程离散化首先要将微分方程离散化.一般采用以下几种方法一般采用以下几种方法:(1)(1)用差商近似导数用差商近似导数 建立数值解法的常用方法建立数值解法的常用方法第8页，本讲稿共69页(2)(2)用数值积分近似积分用数值积分近似积分实际上是矩形法实际上是矩形法宽宽高高第9页，本讲稿共69页(3)(3)用用TaylorTaylor多项式近似并可估计误差多项式近似并可估计误差第10页，本讲稿共69页2 简单的数值方法简单的数值方法 Euler方法的基本原理方法的基本原理将将 在点在点 处进行处进行Taylor展开展开略去略去 项：项：然后用然后用

6、 代替代替 ，即得，即得称上述公式为向前称上述公式为向前Euler 公式。公式。一、一、Euler方法方法第11页，本讲稿共69页若将若将 在点在点 处进行处进行Taylor展开展开略去略去 项：项：然后用然后用 代替代替 ，即得，即得称上述公式为向后称上述公式为向后Euler 公式。公式。向后向后Euler 公式为公式为隐式隐式格式，需要利用格式，需要利用迭代法迭代法求解求解第12页，本讲稿共69页 Euler方法的方法的几何意义几何意义第13页，本讲稿共69页Y=y(x)ab第14页，本讲稿共69页解：解：向前向前Euler公式：公式：例例1：分别利用向前和向后分别利用向前和向后Euler

7、方法方法求解初值问题求解初值问题的的数值数值解解（取步长为（取步长为 ）向后向后Euler公式：公式：第15页，本讲稿共69页具体计算结果具体计算结果:第16页，本讲稿共69页第17页，本讲稿共69页利用数值积分将微分方程离散化利用数值积分将微分方程离散化得梯形公式得梯形公式:解决方法：有的可化为显格式，但有的不行解决方法：有的可化为显格式，但有的不行梯形方法为隐式算法梯形方法为隐式算法二、二、改进的改进的Euler方法方法第18页，本讲稿共69页梯形公式比梯形公式比EulerEuler法精度高一些法精度高一些,但计算量较大但计算量较大 实际计算中只迭代一次，这样建立的预测实际计算中只迭代一次

8、，这样建立的预测校正校正系统称作改进的系统称作改进的Euler公式。公式。第19页，本讲稿共69页第20页，本讲稿共69页例例解解第21页，本讲稿共69页Euler近似解近似解精确解精确解0 1.0.1 1.10.2 1.191820.3 1.277440.4 1.358210.5 1.435130.6 1.50897y0-1y0.1-1.09545y0.2-1.18322y0.3-1.26491y0.4-1.34164y0.5-1.41421y0.6-1.483240 1.0.1 1.097740.2 1.187570.3 1.271290.4 1.350130.5 1.424990.6 1

9、.49657改进改进Euler近似解近似解结果比较结果比较第22页，本讲稿共69页三、三、常微分方程数值解法的常微分方程数值解法的稳定性稳定性设一个数值方法以定步长设一个数值方法以定步长 求解实验方程求解实验方程得到线性差分方程的解得到线性差分方程的解 。当时。当时 ，若，若 ，则称该方法对步长为则称该方法对步长为绝对稳定绝对稳定的；否则称为不稳定的。的；否则称为不稳定的。将数值方法应用于实验方程，若对一切将数值方法应用于实验方程，若对一切都是绝对稳定的，则称区域都是绝对稳定的，则称区域 为该方法的为该方法的绝对稳定域绝对稳定域。上述定义表明，若数值方法可使任何一步产生的误差在后面的上述定义表

10、明，若数值方法可使任何一步产生的误差在后面的计算中都能逐步削弱，则该方法为绝对稳定。计算中都能逐步削弱，则该方法为绝对稳定。第23页，本讲稿共69页例如例如，对于向前，对于向前Euler法：法：将其应用于将其应用于实验实验方程方程当当 时，误差将逐步减弱，故此时时，误差将逐步减弱，故此时方法稳定。方法稳定。向前向前Euler法法绝对稳定域绝对稳定域：当当 因有误差变为因有误差变为 时，时，则有则有第24页，本讲稿共69页四、四、单步方法单步方法的局部误差的局部误差和和阶阶单步法单步法的一般形式的一般形式隐式隐式单步法单步法通常称通常称 为为增量增量函数函数显式显式单步法单步法称称 为某方法在点

11、为某方法在点 的的整体截断整体截断误差误差设设 是是准确准确的，用某种方法计算的，用某种方法计算 时产时产生的截生的截断误差，称为该方法的断误差，称为该方法的局部局部截断截断误差，即误差，即(单步法：在计算单步法：在计算yn+1 时只利用时只利用yn)第25页，本讲稿共69页其中其中 为自然数，则称该方法是为自然数，则称该方法是 阶的或具有阶的或具有 阶精阶精度。度。如果给定方法的如果给定方法的局部截断局部截断误差为误差为如果一个如果一个 阶单步方法的阶单步方法的局部局部截断截断误差为误差为则称则称 为该方法的局部截断误差的为该方法的局部截断误差的主项主项。如向前如向前Euler方法的方法的局

12、部局部截断截断误差误差一阶一阶方法方法第26页，本讲稿共69页 Euler方法的误差分析方法的误差分析对初值问题中的微分方程两端在区间对初值问题中的微分方程两端在区间 上积分上积分如果用如果用左矩形公式左矩形公式计算右端积分，并令计算右端积分，并令其中其中上述等式中如果用上述等式中如果用 代替代替 ，即得向前，即得向前Euler格式。格式。其其局部截断局部截断误差为误差为第27页，本讲稿共69页设设 关于关于 和和 均满足均满足Lipschitz条件，即条件，即和和第28页，本讲稿共69页其中其中而而整体截断整体截断误差为误差为第29页，本讲稿共69页 注意到注意到第30页，本讲稿共69页对于

13、初值问题对于初值问题 ，如果，如果 关于关于 满足满足（向前向前Euler方法的方法的整体截断整体截断误差误差）Lipschitz条件，条件，为对应的为对应的Lipschitz常数，当常数，当时，向前时，向前Euler方法的数值解方法的数值解 一致收敛一致收敛于初值问题于初值问题 的精确解，且的精确解，且整体截断整体截断误差满足估计式误差满足估计式如果如果 ，Euler方法的方法的整体整体截断截断误差为误差为 第31页，本讲稿共69页一、一、Runge-Kutta方法的基本思想方法的基本思想3 龙格龙格-库塔（库塔（Runge-Kutta）方法）方法显式显式单步法的一般形式：单步法的一般形式：

14、R-K方法是利用一些点的线性组合方法是利用一些点的线性组合构造构造增量函数增量函数，使得相应方法的使得相应方法的局部局部截断误差的截断误差的阶数阶数尽可能尽可能高高。二阶二阶Runge-Kutta方法方法确定参数确定参数 ，使得，使得与与 在点在点 的的Taylor展开式有尽可能多的展开式有尽可能多的相同项相同项。第32页，本讲稿共69页比较两式的比较两式的相同项相同项得得方程组有无穷多解方程组有无穷多解第33页，本讲稿共69页若取其一组解若取其一组解则得到则得到改进改进的的Euler公式（公式（二阶二阶方法）方法）若取其另一组解若取其另一组解则得到则得到二阶二阶的的Heun（休恩）公式。（休

15、恩）公式。第34页，本讲稿共69页二、显式二、显式Runge-Kutta方法及其稳定性方法及其稳定性和和设设 是一个正整数，代表使用函数值是一个正整数，代表使用函数值 的个数，的个数，是一些特定是一些特定的权因子（均为的权因子（均为实数），则称下列方法（公式）实数），则称下列方法（公式）为初值问题为初值问题 的的m级级显式显式RungeKutta公式，公式，其中其中第35页，本讲稿共69页类似前面的处理方法，可以得到类似前面的处理方法，可以得到四级四级方法：方法：m=4局部截断局部截断误差误差最常用的一种最常用的一种四阶四阶方法：方法：经典显式经典显式Runge-Kutta公式公式第36页，本

16、讲稿共69页解：解：例例2：用用经典的四阶经典的四阶Runge-Kutta方法求解下列初值问题方法求解下列初值问题 。经典的四阶经典的四阶Runge-Kutta公式：公式：第37页，本讲稿共69页第38页，本讲稿共69页第39页，本讲稿共69页注注：对于显式对于显式N级级R-K方法，最多只能得到方法，最多只能得到N阶方法。阶方法。上述方法的缺陷：计算非常复杂。上述方法的缺陷：计算非常复杂。可通过积分方法确定参数。可通过积分方法确定参数。例例2：确定如下三级三阶显式确定如下三级三阶显式Runge-Kutta公式中的参数：公式中的参数：解：解：对微分方程对微分方程 两边积分得两边积分得第40页，本

17、讲稿共69页采用采用Simpson公式计算上式右端积分项公式计算上式右端积分项可设参数可设参数则有则有选择选择剩余参数剩余参数，使得，使得第41页，本讲稿共69页取取第42页，本讲稿共69页第43页，本讲稿共69页取取利用利用Taylor展开式展开式第44页，本讲稿共69页代入代入当当 时，时，第45页，本讲稿共69页例例3：求求经典四阶的经典四阶的R-K方法的方法的绝对稳定域。绝对稳定域。解：解：第46页，本讲稿共69页其其绝对稳定域绝对稳定域为为三、隐式三、隐式Runge-Kutta方法方法m级级隐式隐式RK方法的一般形式方法的一般形式其中系数的确定方法同其中系数的确定方法同显式显式RK方

18、法完全类似方法完全类似第47页，本讲稿共69页(1)(1)一一级级二阶二阶的隐式的隐式中点中点方法：方法：(2)(2)二二级级四阶四阶的隐式的隐式R-K方法：方法：N级隐式级隐式R-K法法可以达到可以达到2N阶阶缺陷：需要求解缺陷：需要求解非线性方程（组）非线性方程（组）第48页，本讲稿共69页一、一、k步线性步线性多步法多步法4 线性线性多步法多步法/*Linear Mutistep Method and Predictor-Corrector Format*/所谓的线性所谓的线性多步法多步法，指的是某一步解的公式不仅与前，指的是某一步解的公式不仅与前一步的值有关，而且与前面若干步解的值有关

19、的方法。一步的值有关，而且与前面若干步解的值有关的方法。对初值问题对初值问题 两边两边积分积分得得第49页，本讲稿共69页将将 换为节点换为节点取节点取节点 ，构造，构造 的的k+1个点的个点的Lagrange插值多项式：插值多项式：多步多步显式显式公式公式第50页，本讲稿共69页其中其中记记若函数值若函数值 已知，则得已知，则得r+1步步显式显式方法方法第51页，本讲稿共69页如如 时，可得二步时，可得二步显式显式阿达姆斯（阿达姆斯（Adams）格式）格式其中其中第52页，本讲稿共69页 Adams显式公式的显式公式的局部局部截断误差：截断误差：由由Lagrange插值余项知插值余项知其中其

20、中（第（第二二积分积分中值中值定理）定理）k阶方法阶方法第53页，本讲稿共69页取节点取节点 ，构造，构造 的的k+1个点的个点的Lagrange插值多项式：插值多项式：多步多步隐式隐式公式公式第54页，本讲稿共69页其中其中记记则得到则得到r+1步步q+1阶的阶的隐式隐式方法方法如如 时，可得二步时，可得二步隐式隐式阿达姆斯（阿达姆斯（Adams）格式）格式梯形梯形公式公式第55页，本讲稿共69页 常用的一种常用的一种预测预测-校正校正公式：公式：四阶四阶Adams预测预测-校正公式：校正公式：（显式显式）（隐式隐式）初始迭代值由初始迭代值由4阶阶R-K方法计算方法计算第56页，本讲稿共69

21、页例例4：用用Adams预测预测-校正公式校正公式求解下列初值问题求解下列初值问题 。解：解：Adams预测预测-校正公式：校正公式：第57页，本讲稿共69页R-K方法方法Adams预预-校法校法 精确精确解解 0 11.0000000000 0.11.0954461.0954451153 0.21.1832171.1832159566 0.31.2649121.2649110640 0.41.34164135711.3416407864 0.51.41421383341.4142135623 0.61.48323982421.4832396974 0.71.54919338041.54919

22、33384 0.81.61245153641.6124515496 0.91.67331999931.6733200530 1.01.73205071981.7320508075第58页，本讲稿共69页 5 一阶方程组一阶方程组与与高阶高阶方程的数值解法方程的数值解法一、一、一阶微分方程组初值问题的一阶微分方程组初值问题的一般形式一般形式初始条件初始条件:第59页，本讲稿共69页写成写成向量向量的的形式形式:第60页，本讲稿共69页 n=2对应的对应的Runge-Kutta公式公式第61页，本讲稿共69页第62页，本讲稿共69页例例 考虑考虑LorenzLorenz模型：模型：其中参数b=8/

23、3，c=10，r=28解：第63页，本讲稿共69页图中，图中，x x1 1的图形为实线的图形为实线(蓝），蓝），x x2 2的图形为的图形为“*”线（绿），线（绿），x x3 3的图形为的图形为“+”线（红）。线（红）。MatlabMatlab软件计算数值解，计算结果如下图软件计算数值解，计算结果如下图若自变量区间取若自变量区间取0，10，计算结果如下：，计算结果如下：右图的坐标为右图的坐标为x x1 1，x x2 2，x x3 3。左图的坐标为左图的坐标为t t，x x1 1或或x x2 2或或x x3 3。第64页，本讲稿共69页曲线呈震荡发散状三维图形的混沌状若自变量区间取若自变量区间取

24、0，20、0，40，计算结果如下：，计算结果如下：第65页，本讲稿共69页观察结果：观察结果：1、该曲线包含两个“圆盘”，每一个都是由螺线形轨道构成。某些轨道几乎是垂直地离开圆盘中一个而进入另一个。2、随着t的增加，x(t)先绕一个圆盘几圈，然后“跳”到另一个圆盘中。绕第二个圆盘几圈，又跳回原来的圆盘。并以这样的方式继续下去，在每个圆盘上绕的圈数是随机的。第66页，本讲稿共69页作下列作下列变量代换变量代换可将其化为可将其化为一阶方程组一阶方程组的初值问题的初值问题:二、二、高阶高阶微分方程初值问题的微分方程初值问题的一般形式一般形式第67页，本讲稿共69页可类似可类似前述方法前述方法构造构造相应的数值方法求解相应的数值方法求解第68页，本讲稿共69页希望大家考出好成绩!预祝大家研究生涯愉快！第69页，本讲稿共69页