《分析09-常微方程数值解法》由会员分享，可在线阅读，更多相关《分析09-常微方程数值解法（81页珍藏版）》请在金锄头文库上搜索。
1、第八章 常微分方程数值解法,8-1,第八章,常微分方程数值解法,第八章 常微分方程数值解法,8-2,第八章目录,1 欧拉（Euler）方法1.1 Euler法及其简单改进1.2 改进的Euler法 2 龙格库塔（Runge-kutta）方法2.1 龙格-库塔方法的基本思想2.2 二阶龙格-库塔公式2.3 高阶R-K公式2.4 变步长R-K法 3 线性多步法 4 一阶方程组与高阶方程初值问题 5 收敛性与稳定性,第八章 常微分方程数值解法,8-3,第八章 序,许多科学技术问题，例如天文学中的星体运动，空间 技术中的物体飞行，自动控制中的系统分析，力学中的振 动，工程问题中的电路分析等，都可归结为
2、常微分方程的 初值问题。,所谓初值问题，是函数及其必要的导数在积分的起始 点为已知的一类问题，一般形式为：,我们先介绍 简单的一阶问题：,第八章 常微分方程数值解法,8-4,第八章 序,由常微分方程的理论可知：上述问题的解唯一存在。常微分方程求解求什么？应求一满足初值问题（81） 的解函数y = y(x) ，如对下列微分方程：,高等数学中，微分方程求解，如对一阶微分方程： y =f(x,y)是求解解函数y = y(x) ，使满足上述方程。但能够 求出准确的解析函数y(x)的微分方程是很少的，高数 中研究微分方程的求解，是分门别类讨论，对不同类型的 微分方程，求解方法不一样，因此，要求解微分方程
3、，首 先必须认清类型。,第八章 常微分方程数值解法,8-5,微分方程 数值解,而常微分方程 初值问题的数值解法，是要寻求解函数y(x) 在一系列点y(xi ) （离散点）:,上 y(xi )的近似值yi（ i=1,2,n），并且还可由这些（xi,yi) （i=1,2,n）构造插值函数作为近似函数。上述离散点相 邻两点间的距离hi=xi-1-xi 称为步长，若hi 都相等为一定数 h, 则称为定步长，否则为变步长。,由于在实际问题和科学研究中遇到的微分方程往往很 复杂，绝大多数很难，甚至不可能求出解析函数y(x)，因 此只能考虑求其数值解。,本章重点讨论如下 一阶微分方程：,在此基础上介绍一阶微
4、分方程组与 高阶微分方程的数值解法。,第八章 常微分方程数值解法,8-6,1 欧拉（Euler）法,以Euler法及其改进方法为例,说明 常微分方程初值问题数值解法的一般概 念，Euler法很简单，准确度也不高， 介绍此方法的目的，是由于对它的分析 讨论能够比较清楚地显示出方法的一些 特点，而这些特点及基本方法反映了其 它方法的特点。,Euler法用于求 解一阶微分方 程初值问题：,第八章 常微分方程数值解法,8-7,1.1 Euler法及其简单改进,Euler公 式为：,由x0出发x1,x2,xN，而利用此式可算出对应的 y1,y2,yN，式（8-2）称为差分方程（序列yn满足的方 程）。,
5、下面是Euler公式的推导 ：,一、从几何意义出发：y=f (x,y)的解函数y=y (x) 在xoy平 面上是一族解曲线， 而初值问题则是其中一条积分曲线， 假定y = y(x)的曲线如图8-1从给定的点P0(x0,y0) 出发，以P0为切点，作切线，切线斜率为曲线y(x)的切线斜率 y =f (x0,y0)，因此可得切线：（点斜式）,第八章 常微分方程数值解法,8-8,Euler公式的推导(续1）,几何意义：用折线近似解曲线y = y(x)，折线不会偏离太远，因为每项以f (x, y)（斜率）修正。,切线与x = x1交于P1(x1,y1)，在x0,x1上以切线,近似曲线，,第八章 常微分
6、方程数值解法,8-9,Euler公式的推导（续2）,二、利用Taylor级数：将y(x)在xn处展开：,第八章 常微分方程数值解法,8-10,Euler公式的推导（续3）,公式（8-3）称为后退Euler公式,所谓局部载断误差是指以yn代替y(xn)而得到y(xn+1)的近 似值yn+1的误差（只估计这一步的误差）。,三、利用数值微分公式：利用两点公式,后退Euler公式的 局部截断误差为:,第八章 常微分方程数值解法,8-11,Euler公式的推导（续4）,第八章 常微分方程数值解法,8-12,Euler公式的推导（续5）,四、利用数值积分公式：在xn,xn+1上对y(x)=f (x,y(x
7、) 积分,对右端积分项采用不同的数值积分公式，便可得到各种 不同的求解dE初值问题的计算公式。如，以矩形面积代替曲边梯形面积,1）以左矩形面积代替曲边梯形面积如图8-2，亦即以,图8-2,第八章 常微分方程数值解法,8-13,3）以梯形公式（面积）代替曲边形如图8-4则有,式（8-5）称为求dE初值问题的梯形公式，梯形公式看作 是以（xn,yn）(xn+1,yn+1)构造的插值多项式代替被积函数 得到的，而Euler公式则是以左端点函数值近似被积函数 而得到，还可以用多个点做插值多项式近似被积函数构造 另一些精度更高的解微分方程的数值公式，梯形公式比Euler公式更准确一些，误差更小。,Eul
8、er公式的推导（续6）,2）以右矩形面积代替曲边梯形：如图8-3亦即以,第八章 常微分方程数值解法,8-14,Euler公式注释,注1 ：Euler公式为显式，右矩形，梯形公式为隐式；,注2：Euler公式，梯形，右矩形公式为单步法，计算yn+1 只用yn，而中点法公式为多步法（还可如上二所述，构造 多步法）即必须已知yn-1，yn才 能计算yn+1，（求y0 ,y1不能 用此公式。y0 , y1称为多步法的开始值，y0给定， 而y1必须 由其它公式算出，然后才能用中点法）；,注3：前面已有Euler法的局部截断误差：,后退Euler法的局部截断误差：,误差阶：如果局部截断误差,则称方法为P阶
9、的。,第八章 常微分方程数值解法,8-15,显然，步长h越小，阶数P越高，局部截断误差越小，当 然计算精度越高；,注4：梯形法是几阶？梯形法精度比Euler法高，阶数肯定 比Euler法高，其实我们可以利用数值积分公式的误差估 计式，因为我们是用梯形数值积 分公式计算,因此由积分中梯形公式的误差知此,时的局部截断误差为：,梯形法为2阶方法 ！,Euler法，后退Euler法为1阶方法，而中点法为2 阶，,Euler公式注释（续）,第八章 常微分方程数值解法,8-16,关于Euler法的整体截断误差注释,注5：关于Euler法 的整体截断误差：,Euler方法的局部截断误差公式为：,实际计算时，
10、yn是y(xn) 的近似值，因此，计算过程 中除每步所产生的局部截断 误差外，还有因前面的计算不准确而引起的误差。在不考 虑舍入误差的情况下，称y(xn+1)与yn+1之差为整体截断误差 ，记为：,下面讨论Euler方法的整体截断误差。,为简便起见，假定函数f (x, y)充分光滑，问题（8-1）解 y(x)在a, b上二阶连续可微，于是由式（8-6），局部截断 误差有界，即存在M0， 使得对任意xa,b，都有 |y(x)|M，从而有：,（紧接下屏）,第八章 常微分方程数值解法,8-17,式（8-7）表明，Euler方法的整体截断误差与h同阶，当h0时，en0。,关于Euler法的整体截断误差
11、注释（续）,反复递推得：,第八章 常微分方程数值解法,8-18,结 论,对于实际问题来说，由于L，M 难以估计，因此（8-6） 很难应用，而且上述推导过程中一再放大了误差上限，这 样的估计往往也很保守，远远大于实际的误差，但是，从 估计式（8-7）却可以得到下面很有用处的结论。,1）当h0时，en0即，,亦即数值解yn，一致收敛于初值问题（8-1）的真解y(xn)，,并且，Euler法的整体截断误差的阶为O(h)与h同阶，比局 部截断误差低一阶。,2）舍入误差局部截断误差 对实际计算结果有影响，并且随h减少 而减少或增大。,第八章 常微分方程数值解法,8-19,3）计算结果与解法的阶数p，真解
12、的导数y (p+1)有关，p越大，h p+1越小， |y(p+1)()|的上限越大，M 也越大，因此为保证精度当然应选阶数p较高的方法。但如果M 很大，当f (x,y)是分段连续的函数时，则应采用低阶的方法如用Euler法。,结 论（续）,4）计算结果还与开始值的精度有关，为使这种误差的影响不致于超过局部 截断误差，对多步法，应采用跟多步法同阶的方法计算开始值。,第八章 常微分方程数值解法,8-20,1.2 改进的Euler法,梯形公式为二阶方法，但却是隐式格式，即若利用梯 形公式求yn+1，就要求解方程（8-5）式，计算量较大，通 常在实际计算时，将Euler法与梯形公式合起来使用，即 先使
13、用Euler公式,由（xn,yn）算出yn+1，记为yn+1(0)，称为 预测值，然后用梯形公式去提高精度，将yn+1(0) 校正为较 准确的值：,由于函数f (x,y)满足Lipschitz条件，容易得出:,第八章 常微分方程数值解法,8-21,改进的Euler法(续）,第八章 常微分方程数值解法,8-22,预测校正型公式,实际经验表明，式（8-8）的迭代效果主要体现在第一次， 由此构成如下的预测校正型公式:,此式称为改进的Euler公式，为上机计算编程方便，常将式 （8-9）改写为:,下面分析改进的Euler 公式的局部截断误差：,第八章 常微分方程数值解法,8-23,改进的Euler公式
14、的局部截断误差分析,假定yn = y(xn) ，y(xn+1) 的 Taylor展式为:,对于改进的Euler公式，由于,这说明改进的Euler法的局部截断误差为O(h3)，比Euler公式高一阶，是二阶方法。,第八章 常微分方程数值解法,8-24,改进的Euler公式举例,例1,这些结果在表8-1中，可见计算结果的精度，Euler法与 后退Euler法差不多，与准确值相比较Euler法偏小，而后 退Euler法偏大；中点法与梯形法精度同为2阶，但梯形法 更好一些，这跟它们局部截断误差的符号，阶数和系 数的大小是完全一致的。,表见下屏：,第八章 常微分方程数值解法,8-25,表格8-1,表8-
15、1 y = y, y(0)=1 的数值解(h=0.1),第八章 常微分方程数值解法,8-26,表格8-2,而表8-2是分别取了不同的h=0.1 ,h=0.01，h=0.001, h=0.0001，还是利用这些公式，经过若干步的计算（h越 小，计算量越大）算到y(1)的近似值，可见：随着h的减小，y(1)的近似值的精度在提高，0.01比 0.001差，即0.001比0.01时的y(1)准确。,（紧接下屏）,第八章 常微分方程数值解法,8-27,表8-2计算结果说明（续）,但h太小，到h=0.0001时却又变得误差大了，这与前面所说h越小，p阶越高，应该局部截断误差越小，因而计算精度更高矛盾了，为
16、什么会产生这种情况呢？这是由于h太小而引起计算量大因而造成了舍入误差和截断误差的积累，这种情况由于初值问题不同可能会影响更大，偏离更严重，如下面的例2 。这种问题实际上是稳定性问题，我们将会讨论方法的稳定性，由此得出对h有一定的要求的稳定性制区域。,第八章 常微分方程数值解法,8-28,例2,求解初值问题y=-20y，y(0)=1，计算y(1)的近似值。,解：类似于例1，用欧拉法、后退欧拉法、中点法、梯形法求解，得到如下表8.3。 表8.3 y = 20,y(0)= 1的解y(1)的近似值(y(1)=0.206116E8),由表8.3可见，尽管中点法的阶数与梯形法相同，比 欧拉法和后退欧拉法的阶数高，计算结果的精度却很糟 糕。此外，尽管欧拉法与后退欧拉法的阶数相同，但欧 拉法计算结果的精度，当h=0.1时却比后退欧拉法差。,
