第九章常微方程数值解法

第九章常微方程数值解法第章9章常微分方程数值解法8-1第1页，课件共31页，创作于2023年2月第9章目录§1欧拉（Euler）方法

1.1Euler法及其简单改进1.2改进的Euler法§2龙格库塔（Runge-kutta）方法

2.1龙格-库塔方法的基本思想2.2二阶龙格-库塔公式2.3高阶R-K公式2.4变步长R-K法§3线性多步法§4一阶方程组与高阶方程初值问题§5收敛性与稳定性第2页，课件共31页，创作于2023年2月第8章序

许多科学技术问题，例如天文学中的星体运动，空间技术中的物体飞行，自动控制中的系统分析，力学中的振动，工程问题中的电路分析等，都可归结为常微分方程的初值问题。

所谓初值问题，是函数及其必要的导数在积分的起始点为已知的一类问题，一般形式为：

我们先介绍简单的一阶问题：第3页，课件共31页，创作于2023年2月第9章序由常微分方程的理论可知：上述问题的解唯一存在。

常微分方程求解求什么？应求一满足初值问题（8—1）的解函数y=y(x)，如对下列微分方程：《高等数学》中，微分方程求解，如对一阶微分方程：y

=f(x,y)是求解解函数y=y(x)，使满足上述方程。但能够求出准确的解析函数y(x)的微分方程是很少的，《高数》中研究微分方程的求解，是分门别类讨论，对不同类型的微分方程，求解方法不一样，因此，要求解微分方程，首先必须认清类型。第4页，课件共31页，创作于2023年2月微分方程数值解而常微分方程初值问题的数值解法，是要寻求解函数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)，因此只能考虑求其数值解。本章重点讨论如下一阶微分方程：在此基础上介绍一阶微分方程组与高阶微分方程的数值解法。第5页，课件共31页，创作于2023年2月§1欧拉（Euler）法

以Euler法及其改进方法为例,说明常微分方程初值问题数值解法的一般概念，Euler法很简单，准确度也不高，介绍此方法的目的，是由于对它的分析讨论能够比较清楚地显示出方法的一些特点，而这些特点及基本方法反映了其它方法的特点。Euler法用于求解一阶微分方程初值问题：第6页，课件共31页，创作于2023年2月1.1Euler法及其简单改进

Euler公式为：由x0出发

x1,x2,…,xN，而利用此式可算出对应的y1,y2,…,yN，式（8-2）称为差分方程（序列{yn}满足的方程）。下面是Euler公式的推导：一、从几何意义出发：y

=f(x,y)的解函数y=y(x)在xoy平面上是一族解曲线，而初值问题则是其中一条积分曲线假定y=y(x)的曲线如图8-1从给定的点P0(x0,y0)出发，以P0为切点，作切线，切线斜率为曲线y(x)的切线斜率

y

=f(x0,y0)，因此可得切线：（点斜式）P1P2y(x)P0

x2x1x0第7页，课件共31页，创作于2023年2月Euler公式的推导(续1）几何意义：用折线近似解曲线y=y(x)，折线不会偏离太远，因为每项以f(x,y)（斜率）修正。切线与x=x1交于P1(x1,y1)，在[x0,x1]上以切线近似曲线，

第8页，课件共31页，创作于2023年2月Euler公式的推导（续2）二、利用Taylor级数：将y(x)在xn处展开：

第9页，课件共31页，创作于2023年2月Euler公式的推导（续4）第10页，课件共31页，创作于2023年2月Euler公式的推导（续5）四、利用数值积分公式：在[xn,xn+1]上对y

(x)=f(x,y(x))积分

对右端积分项采用不同的数值积分公式，便可得到各种不同的求解dE初值问题的计算公式。如，以矩形面积代替曲边梯形面积1）以左矩形面积代替曲边梯形面积如图8-2，亦即以yf(x,y)xnxxn+1图8-2

第11页，课件共31页，创作于2023年2月yf(x,y)xnxxn+1图8-3yf(x,y)xnxxn+1图8-43）以梯形公式（面积）代替曲边形如图8-4则有

式（8-5）称为求dE初值问题的梯形公式，梯形公式看作是以（xn,yn）(xn+1,yn+1)构造的插值多项式代替被积函数得到的，而Euler公式则是以左端点函数值近似被积函数而得到，还可以用多个点做插值多项式近似被积函数构造另一些精度更高的解微分方程的数值公式，梯形公式比Euler公式更准确一些，误差更小。Euler公式的推导（续6）2）以右矩形面积代替曲边梯形(后退的梯形公式）：如图8-3

第12页，课件共31页，创作于2023年2月Euler公式注释注1：Euler公式为显式，后退的Euler公式，梯形公式为隐式；注2：Euler公式，梯形，后退的Euler公式为单步法，计算yn+1只用yn，而中点法公式为多步法（还可如上二所述，构造多步法）即必须已知yn-1，yn才能计算yn+1，（求y0,y1不能用此公式。y0,y1称为多步法的开始值，y0给定，而y1必须由其它公式算出，然后才能用中点法）；注3：前面已有Euler法的局部截断误差：后退Euler法的局部截断误差：

误差阶：如果局部截断误差则称方法为P阶的。第13页，课件共31页，创作于2023年2月显然，步长h越小，阶数P越高，局部截断误差越小，当然计算精度越高；

注4：梯形法是几阶？梯形法精度比Euler法高，阶数肯定比Euler法高，其实我们可以利用数值积分公式的误差估计式，因为我们是用梯形数值积分公式计算因此由积分中梯形公式的误差知此时的局部截断误差为：∴梯形法为2阶方法！Euler法，后退Euler法为1阶方法，而中点法为2阶，Euler公式注释（续）第14页，课件共31页，创作于2023年2月关于Euler法的整体截断误差注释注5：关于Euler法的整体截断误差：

实际计算时，yn是y(xn)的近似值，因此，计算过程中除每步所产生的局部截断误差外，还有因前面的计算不准确而引起的误差。在不考虑舍入误差的情况下，称y(xn+1)与yn+1之差为整体截断误差，记为：下面讨论Euler方法的整体截断误差。

为简便起见，假定函数f(x,y)充分光滑，问题（8-1）解y(x)在[a,b]上二阶连续可微，于是由式（8-6），局部截断误差有界，即存在M>0，使得对任意x

[a,b]，都有|y′(x)|≤M，从而有：（紧接下屏）第15页，课件共31页，创作于2023年2月结论

对于实际问题来说，由于L，M难以估计，因很难应用，而且上述推导过程中一再放大了误差上限，这样的估计往往也很保守，远远大于实际的误差，但是，从估计式（8-7）却可以得到下面很有用处的结论。

1）当h

0时，en

0即，亦即数值解yn，一致收敛于初值问题（8-1）的真y(xn)，并且，Euler法的整体截断误差的阶为O(h)与h同阶，比局部截断误差低一阶。2）舍入误差局部截断误差对实际计算结果有影响，并且随h减少而减少或增大。第16页，课件共31页，创作于2023年2月3）计算结果与解法的阶数p，真解的导数y(p+1)有关，p越大，hp+1越小，|y(p+1)(ξ)|的上限越大，M也越大，因此为保证精度当然应选阶数p较高的方法。但如果M很大，当f(x,y)是分段连续的函数时，则应采用低阶的方法如用Euler法。结论（续）

4）计算结果还与开始值的精度有关，为使这种误差的影响不致于超过局部截断误差，对多步法，应采用跟多步法同阶的方法计算开始值。第17页，课件共31页，创作于2023年2月1.2改进的Euler法

梯形公式为二阶方法，但却是隐式格式，即若利用梯形公式求yn+1，就要求解方程（8-5）式，计算量较大，通常在实际计算时，将Euler法与梯形公式合起来使用，即先使用Euler公式,由（xn,yn）算出yn+1，记为yn+1(0)，称为

预测值，然后用梯形公式去提高精度，将yn+1(0)

校正为较准确的值：由于函数f(x,y)满足Lipschitz条件，容易得出:第18页，课件共31页，创作于2023年2月改进的Euler法(续）第19页，课件共31页，创作于2023年2月预测──校正型公式

实际经验表明，式（8-8）的迭代效果主要体现在第一次，由此构成如下的预测──校正型公式:此式称为改进的Euler公式，为上机计算编程方便，常将式

（8-9）改写为:下面分析改进的Euler公式的局部截断误差：

第20页，课件共31页，创作于2023年2月改进的Euler公式的局部截断误差分析假定yn=y(xn)，y(xn+1)的Taylor展式为:对于改进的Euler公式，由于

这说明改进的Euler法的局部截断误差为O(h3)，比Euler公式高一阶，是二阶方法。第21页，课件共31页，创作于2023年2月改进的Euler公式举例例1

这些结果在表8-1中，可见计算结果的精度，Euler法与后退Euler法差不多，与准确值相比较Euler法偏小，而后退Euler法偏大；中点法与梯形法精度同为2阶，但梯形法更好一些，这跟它们局部截断误差的符号，阶数和系

数的大小是完全一致的。表见下屏：第22页，课件共31页，创作于2023年2月表格8-1表8-1

y

=

y,y(0)=1的数值解(h=0.1)

x精确解欧拉法后退欧拉中点法梯形法.1.904837.900000.909091.900000.904762.2.808731.810000.826446.820000.818594.3.740818.729000.751315.736000.740633.4.670320.656100.683013.627800.670096.5.060531.590490.620921.601440.606278.6.548812.531441.654474.552512.548537.7.496585.478298.5131458.490938.496295.8.449329.430467.466507.454324.449029.9.406570.387421.424098.400073.4062641.367879.348679.385543.374310.367573第23页，课件共31页，创作于2023年2月表格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)准确。Y′=-y,y(0)=1的解y(1)的近似值(y(1)=0.367879)h欧拉法后退欧拉法中点法梯形法0.1.348678.385543.374310.3675730.01.366033.369711.367944.3678770.001.367700.368052.367879.3678760.0001.367800.367800.367881.368020（紧接下屏）第24页，课件共31页，创作于2023年2月表8-2计算结果说明（续）

但h太小，到h=0.0001时却又变得误差大了，这与前面所说h越小，p阶越高，应该局部截断误差越小，因而计算精度更高矛盾了，为什么会产生这种情况呢？这是由于h太小而引起计算量大因而造成了舍入误差和截断误差的积累，这种情况由于初值问题不同可能会影响更大，偏离更严重，如下面的例2。这种问题实际上是稳定性问题，我们将会讨论方法的稳定性，由此得出对h有一定的要求的稳定性制区域。第25页，课件共31页，创作于2023年2月§5收敛性与稳定性

通过前面的讨论可以看到，微分方

程数值解法的基本思想是：

通过某种离散化手段，将微分方程

转化为差分格式求解。从理论上说，这

里需要解决两个问题，一是当步长h

0

时，差分解，这就是差分格式的收敛性

问题；二是利用差分格式求解时，初始

误差及计算过程中的舍入误差能否得到

控制，这就是差分格式的稳定性问题。

本节对这两个问题作一简单介绍。第26页，课件共31页，创作于2023年2月5.1收敛性

用某种差分格式解微分方程时，若对求解区间[a,b]中的任一点x，当h

0时差分解yn