第十五章常微分方程的解法文档格式.docx

1、代替代入（1）中的微分方程，则得化简得如果用的近似值代入上式右端，所得结果作为的近似值，记为，则有 （2）这样，问题（1）的近似解可通过求解下述问题 （3）得到，按式（3）由初值可逐次算出式（3）是个离散化的问题，称为差分方程初值问题。需要说明的是，用不同的差商近似导数，将得到不同的计算公式。）用数值积分方法将问题（1）的解表成积分形式，用数值积分方法离散化。例如，对微分方程两端积分，得 （4）右边的积分用矩形公式或梯形公式计算。）Taylor多项式近似将函数在处展开，取一次Taylor多项式近似，则得再将，得到离散化的计算公式 以上三种方法都是将微分方程离散化的常用方法，每一类方法又可导出不

2、同形式的计算公式。其中的Taylor展开法，不仅可以得到求数值解的公式，而且容易估计截断误差。2 欧拉（Euler）方法2.1 Euler方法Euler 方法就是用差分方程初值问题（3）的解来近似微分方程初值问题（1）的解，即由公式（3）依次算出这组公式求问题（1）的数值解称为向前Euler公式。如果在微分方程离散化时，用向后差商代替导数，即，则得计算公式 （5）用这组公式求问题（1）的数值解称为向后Euler公式。向后Euler法与Euler法形式上相似，但实际计算时却复杂得多。向前Euler公式是显式的，可直接求解。向后Euler公式的右端含有，因此是隐式公式，一般要用迭代法求解，迭代公式

3、通常为 （6）2.2 Euler方法的误差估计 对于向前Euler公式（3）我们看到，当时公式右端的都是近似的，所以用它计算的会有累积误差，分析累积误差比较复杂，这里先讨论比较简单的所谓局部截断误差。假定用（3）式时右端的没有误差，即，那么由此算出 （7）局部截断误差指的是，按（7）式计算由这一步的计算值与精确值之差为了估计它，由Taylor展开得到的精确值是 （8）（7）、（8）两式相减（注意到）得 （9）即局部截断误差是阶的，而数值算法的精度定义为：若一种算法的局部截断误差为，则称该算法具有阶精度。显然越大，方法的精度越高。式（9）说明，向前Euler方法是一阶方法，因此它的精度不高。3

4、改进的Euler方法3.1 梯形公式利用数值积分方法将微分方程离散化时，若用梯形公式计算式（4）中之右端积分，即并用这就是求解初值问题（1）的梯形公式。直观上容易看出，用梯形公式计算数值积分要比矩形公式好。梯形公式为二阶方法。梯形公式也是隐式格式，一般需用迭代法求解，迭代公式为 （10） 由于函数关于满足Lipschitz条件，容易看出其中为Lipschitz常数。因此，当时，迭代收敛。但这样做计算量较大。如果实际计算时精度要求不太高，用公式（10）求解时，每步可以只迭代一次，由此导出一种新的方法改进Euler法。3.2 改进Euler法按式（5）计算问题（1）的数值解时，如果每步只迭代一次，

5、相当于将Euler公式与梯形公式结合使用：先用Euler公式求的一个初步近似值，称为预测值，然后用梯形公式校正求得近似值，即 （11）式（11）称为由Euler公式和梯形公式得到的预测校正系统，也叫改进Euler法。为便于编制程序上机，式（11）常改写成 （12）改进Euler法是二阶方法。4 龙格库塔（RungeKutta）方法回到Euler方法的基本思想用差商代替导数上来。实际上，按照微分中值定理应有注意到方程就有 （13）不妨记，称为区间上的平均斜率。可见给出一种斜率，（13）式就对应地导出一种算法。向前Euler公式简单地取为，精度自然很低。改进的Euler公式可理解为取，的平均值，其

6、中，这种处理提高了精度。如上分析启示我们，在区间内多取几个点，将它们的斜率加权平均作为，就有可能构造出精度更高的计算公式。这就是龙格库塔方法的基本思想。4.1 首先不妨在区间内仍取2个点，仿照（13）式用以下形式试一下 （14）为待定系数，看看如何确定它们使（14）式的精度尽量高。为此我们分析局部截断误差，因为，所以（14）可以化为 （15）在点作了Taylor展开。（15）式又可表为注意到中，可见为使误差，只须令 （16）待定系数满足（16）的（15）式称为2阶龙格库塔公式。由于（16）式有4个未知数而只有3个方程，所以解不唯一。不难发现，若令 ，即为改进的Euler公式。可以证明，在内只取

7、2点的龙格库塔公式精度最高为2阶。4.2 4阶龙格库塔公式要进一步提高精度，必须取更多的点，如取4点构造如下形式的公式： （17）其中待定系数共13个，经过与推导2阶龙格库塔公式类似、但更复杂的计算，得到使局部误差的11个方程。取既满足这些方程、又较简单的一组，可得 （18）这就是常用的4阶龙格库塔方法（简称RK方法）。5 线性多步法以上所介绍的各种数值解法都是单步法，这是因为它们在计算时，都只用到前一步的值，单步法的一般形式是 （19）称为增量函数，例如Euler方法的增量函数为，改进Euler法的增量函数为如何通过较多地利用前面的已知信息，如，来构造高精度的算法计算，这就是多步法的基本思想

8、。经常使用的是线性多步法。让我们试验一下，即利用计算的情况。从用数值积分方法离散化方程的（4）式出发，记，式中被积函数用二节点的插值公式得到（因，所以是外插。 （20）此式在区间上积分可得于是得到 （21）注意到插值公式（20）的误差项含因子，在区间上积分后得出，故公式（21）的局部截断误差为，精度比向前Euler公式提高1阶。若取可以用类似的方法推导公式，如对于有 （22）其局部截断误差为如果将上面代替被积函数用的插值公式由外插改为内插，可进一步减小误差。内插法用的是，取时得到的是梯形公式，取时可得 （23）与（22）式相比，虽然其局部截断误差仍为，但因它的各项系数（绝对值）大为减小，误差还

9、是小了。当然，（23）式右端的未知，需要如同向后Euler公式一样，用迭代或校正的办法处理。6 一阶微分方程组与高阶微分方程的数值解法6.1 一阶微分方程组的数值解法设有一阶微分方程组的初值问题 （24）若记，则初值问题（24）可写成如下向量形式 （25）如果向量函数在区域：满足Lipschitz条件，即存在，使得对，都有那么问题（25）在上存在唯一解问题（25）与（1）形式上完全相同，故对初值问题（1）所建立的各种数值解法可全部用于求解问题（25）。6.2 高阶微分方程的数值解法高阶微分方程的初值问题可以通过变量代换化为一阶微分方程组初值问题。设有阶常微分方程初值问题 （26）引入新变量，问

10、题（26）就化为一阶微分方程初值问题 （27）然后用6.1中的数值方法求解问题（27），就可以得到问题（26）的数值解。最后需要指出的是，在化学工程及自动控制等领域中，所涉及的常微分方程组初值问题常常是所谓的“刚性”问题。具体地说，对一阶线性微分方程组 （28）阶方阵。若矩阵的特征值满足关系则称方程组（28）为刚性方程组或Stiff方程组，称数为刚性比。对刚性方程组，用前面所介绍的方法求解，都会遇到本质上的困难，这是由数值方法本身的稳定性限制所决定的。理论上的分析表明，求解刚性问题所选用的数值方法最好是对步长不作任何限制。7 Matlab解法7.1 Matlab数值解7.1.1 非刚性常微分方

11、程的解法Matlab的工具箱提供了几个解非刚性常微分方程的功能函数，如ode45，ode23，ode113，其中ode45采用四五阶RK方法，是解非刚性常微分方程的首选方法，ode23采用二三阶RK方法，ode113采用的是多步法，效率一般比ode45高。Matlab的工具箱中没有Euler方法的功能函数。）对简单的一阶方程的初值问题改进的Euler公式为我们自己编写改进的Euler 方法函数eulerpro.m如下：function x,y=eulerpro（fun,x0,xfinal,y0,n）;if nargin5,n=50;endh=（xfinal-x0）/n;x（1）=x0;y（1）

12、=y0;for i=1:n x（i+1）=x（i）+h; y1=y（i）+h*feval（fun,x（i）,y（i）; y2=y（i）+h*feval（fun,x（i+1）,y1）; y（i+1）=（y1+y2）/2;例1 用改进的Euler方法求解解 编写函数文件doty.m如下：function f=doty（x,y）;f=-2*y+2*x2+2*x;在Matlab命令窗口输入：x,y=eulerpro（doty,0,0.5,1,10）即可求得数值解。）ode23，ode45，ode113的使用Matlab的函数形式是t,y=solver（F,tspan，y0）这里solver为ode45，ode23，ode113，输入