中北大学《数值分析》常微分方程初值问题的数值解法专业：班级：学号：姓名:日期: 2012.12.26常微分方程初值问题的数值解法摘 要微分方程的数值解法在科学技术及生产实践等多方面应用广泛. 文章分析了构造常微分方程初值问题数值解法的三种常用基本方法，差商代替导数法，数值积分法及待定系数法，推导出了Euler 系列公式及三阶龙格-库塔公式，指出了各公式的优劣性及适用条件，并对Euler 公式的收敛性、稳定性进行了分析。

AbstractThe numerical solution of differential equations is widely used in science, technology, production practices and many other fields. This paper analyzed three kinds of basic methods for constructing numerical solutions for initial value problem of ordinary differential equations ：difference quotient instead of derivative method, numerical integral method and undetermined coefficients method. At the same time, the paper deduces the Euler series formula and the classical third order Runge-Kutta formula. In addition, the paper pointed out the advantages and disadvantages of each formula and application condition, it also analyzed the convergence and stability of the Euler formula.1.引言科学技术及实际生产实践中的许多问题都可归结为微分方程的求解问题，使用较多的是常微分方程初值问题的求解。

对于一阶常微分方程的初值问题000dy /dx f (x,y),y(x )y ,x x b ==<<，其中f 为已知函数，0y 是初始值。

如果函数f 关于变量y 满足Lipschitz 条件，则初值问题有唯一解。

只有当f 是一些特殊类型的函数时，才能求出问题的解析解，但一般情况下都满足不了生产实践与科学技术发展的需要，因此通常求其数值解法。

2.主要算法数值解法是一种离散化的方法，可以求出函数的精确解在自变量一系列离散点处的近似值。

基本思想是离散化，首先要将连续区间离散化，对连续区域[]0x ,b 进行剖分01n 1n x x x x b -<<Λ<<=，n n 1n h x x +=-为步长；其次将其函离散化，基本方法有差商代替导数法、积分插值方法、待定系数法；最后研究其稳定性、收敛性质。

2.1差商法分别用一阶向前、向后、中心差商近似代替y(x)在n y(x)x =处的导数n y'(x )，可以求出Euler 系列公式.2.1.1显式Euler 公式用一阶向前差商代替导数，即n y'(x )≈n 1n y(x )y(x )h+-则n 1n n n n n y(x )y(x )hy '(x )y(x )hf (x ,y(x ))+=+=+因为i i y y(x )≈，公式可以简记为n 1n n n y y hf (x ,y ),n 0,1...+=+=即为显式Euler 公式。

Euler 公式是最简单的一种数值解法，一阶的，精度较差，可直接求解，有明显的几何意义，也称为折线法，但其方法对于更复杂的情况有着较为普遍的意义。

2.1.2隐式Euler 公式用一阶向后差商代替导数，即n n 1n y(x )y(x )y'(x )h+-≈则n n 1n n 1n n y(x )y(x )hy'(x )y(x )hf (x ,y(x ))++=+=+简记为n n 1n n y y hf (x ,y ),n 0,1...-=+=即为隐式Euler 公式或后退Euler 公式，计算比显式麻烦，但稳定性好。

2.1.3两步Euler 公式二阶中心差商代替导数n 1n 1n y(x )y(x )y'(x )2h+--≈则n 1n 1n n y y 2hf (x ,y ),n 0,1...,+-=+=即为两步Euler 公式，稳定性较单步好。

2.2积分差值法 2.2.1梯形公式dy /dx f (x,y)=在区间]n n 1x ,x +⎡⎣上求积分，n 1nx n 1n x y(x )y(x )f (x,y(x))dx ++-=⎰，对右端积分采用积分插值公式n 1nx n n n 1,n 1x hf (x,y(x))dx f (x ,y(x ))f (x y )2+++⎡⎤≈+⎦⎣⎰，即得n 1n n n n 1n 1hy y [f (x ,y )f (x ,y )]2+++=++，为梯形公式，也可看成显式Euler 公式与隐式Euler 公式的算术平均，可以证明梯形方法是二阶方法。

实际计算时，初始近似n 1y +，可由Euler 法求解，则梯形法每步完整的计算公式为(0)n n n n 1(m 1)(m)n n n 1n 1n 1y y hf (x ,y )h y y [f f (x ,y )],m 0,1...2+++++=+=++=⎧⎨⎩又称为改进的Euler 公式。

由Euler 公式给出预测值，再由梯形进行校正，收敛速度较快，一般只需求出两次迭代即可满足要求。

2.2.2中点法同样，对dy /dx f (x,y)=在区间]n 1n 1x ,x -+⎡⎣上求积分，n 1nx n 1n x y(x )y(x )f (x,y(x))dx ++-=⎰并对右端采用Gauss 公式，得n 1n 1n n y y 2hf (x ,y ),n 0,1...+-=+=也称为两步Euler 公式。

对上式右端积分采用Simpson 公式n 1n 1x n 1,n 1n n n 1n 1x hf (x,y(x))dx [f (x y )4f (x ,y )f (x ,y )]3+---++≈++⎰即n 1n 1n 1,n 1n n n 1n 1hy y [f (x y )4f (x ,y )f (x ,y )]3+---++=+++称为隐式两步法。

2.3待定系数法 2.3.1三阶龙格库塔公式R-K 基本思想即用位于]n n 1x ,x +⎡⎣上的若干个点处的线性组合来近似它，将三阶R-K 公式一般形式表示为:10101Y y ,F f (x ,Y )==202112022Y y ha F ,F f (x hc ,Y )=+=+303113223033Y y h(a F a F ),F f (x hc ,Y )=++=+213132123123a ,a ,a ,b ,b ,b ,c ,c ,c 是未知常数，我们的主要目的要目的是估算12F ,F 和F 3，它们分别是00203y'(x ),y'(x hc ),y'(x hc )++的近似解.对导数进行近似，可以从下列正交逼近中得出0010202303y(x h)y(x )h(b y'(x )b y'(x hc )b y'(x hc ))+≈+++++此正交公式可以使多项式的精度达到2阶，待定系数得出如下条件:123b b b 1++= 22331b c b c 2+=2322331b c b c 3+=33221b a c 6=满足四个条件，令21c 2=，得32131321231213c 1,a ,a 0,a ,b ,b 0,b 3344=======即为三阶Kutta 公式。

R-K 公式的推导是基于Taylor 展开，对光滑度要求较高，如果解充分光滑，则结果精度较高。

3.数值实验3.1实验内容科学计算中经常遇到微分方程（组）初值问题，需要利用Euler 法，改进Euler 法，Rung-Kutta 方法求其数值解，诸如以下问题：（1）402(0)1⎧'=-⎪<≤⎨⎪=⎩xy xy y x y分别取h=0.1,0.2,0.4时数值解。

初值问题的精确解=y 。

（2） 2210(1)0y x y x y '⎧=--≤≤⎨-=⎩取步长h=0.1，用四阶标准R-K 方法求值。

3.2实验前的预备知识1、 熟悉各种初值问题的算法；2、 明确各种算法的精度与所选步长有密切关系；3、 通过计算更加了解各种算法的优越性。

3.3实验方法或步骤1、 根据初值问题数值算法，编程计算；2、 试分别取不同步长，考察某节点j x 处数值解的误差变化情况；3、 试用不同算法求解某初值问题，结果有何异常；4、 分析各个算法的优缺点。

3.4实验过程1、402(0)1xy xyy xy⎧'=-⎪<≤⎨⎪=⎩分别取h=0.1,0.2,0.4时数值解。

（1）、步长为0.1时的运行结果：[x,y]=euler('doty',0,2,1,20)[x,y]=euler('doty',0,2,1,20)x =Columns 1 through 130 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000 0.9000 1.0000 1.1000 1.2000Columns 14 through 211.3000 1.4000 1.5000 1.6000 1.7000 1.8000 1.90002.0000y =Columns 1 through 131.0000 1.0000 1.0300 1.0871 1.1648 1.2556 1.3521 1.4485 1.5404 1.6249 1.7002 1.7655 1.8205Columns 14 through 211.8657 1.9019 1.9301 1.9514 1.9672 1.9784 1.9862 1.9915（2）、步长为0.2时的结果：[x,y]=euler('doty',0,2,1,10)[x,y]=euler('doty',0,2,1,10)x =0 0.2000 0.4000 0.6000 0.8000 1.0000 1.2000 1.4000 1.6000 1.8000 2.0000y = 1.0000 1.0000 1.1200 1.3161 1.5229 1.6995 1.8303 1.9155 1.9639 1.9872 1.9964（3）步长为0.4的结果：[x,y]=euler('doty',0,2,1,5)x =0 0.4000 0.8000 1.2000 1.6000 2.0000y =1.0000 1.0000 1.4800 1.8713 1.9991 2.0003改进的Euler法的求解：（1）、步长为0.1时的运行结果：[x,y]=eulerpro('doty',0,2,1,20)x =Columns 1 through 130 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000 0.9000 1.0000 1.1000 1.2000Columns 14 through 211.3000 1.4000 1.5000 1.6000 1.7000 1.8000 1.90002.0000y = Columns 1 through 131.0000 1.0150 1.0575 1.1220 1.2017 1.2898 1.3806 1.4697 1.5540 1.6312 1.7002 1.7603 1.8116Columns 14 through 211.8544 1.8894 1.9174 1.9393 1.9562 1.9690 1.9784 1.9852（2）、步长为0.2时，输入[x,y]=eulerpro('doty',0,2,1,10)[x,y]=eulerpro('doty',0,2,1,10)x =0 0.2000 0.4000 0.6000 0.8000 1.00001.2000 1.4000 1.6000 1.80002.0000y =1.0000 1.0600 1.2045 1.3814 1.5519 1.69551.8050 1.8820 1.9322 1.9629 1.9804（3）步长为0.4的结果：[x,y]=eulerpro('doty',0,2,1,5)x =0 0.4000 0.8000 1.2000 1.6000 2.0000y =1.0000 1.2400 1.5605 1.7797 1.8892 1.9343四阶RK方法：（1）步长为0.1时的运行结果：[x,y]=RK(0,0.1,2,1)x =Columns 1 through 130 0.1000 0.2000 0.3000 0.4000 0.50000.6000 0.7000 0.8000 0.9000 1.0000 1.10001.2000Columns 14 through 211.3000 1.4000 1.5000 1.6000 1.7000 1.80001.90002.0000y =Columns 1 through 131.0000 0.9094 0.8358 0.7772 0.7327 0.70180.6842 0.6802 0.6897 0.7130 0.7500 0.80040.8637Columns 14 through 210.9389 1.0249 1.1203 1.2233 1.3324 1.4571.5619 1.6795（2）步长为0.2的结果[x,y]=RK(0,0.2,2,1)x =0 0.2000 0.4000 0.6000 0.8000 1.0000 1.2000 1.4000 1.6000 1.8000 2.0000y =1.0000 0.8358 0.7327 0.6842 0.6898 0.7500 0.8637 1.0250 1.2234 1.4458 1.6795（3）：步长为0.4的结果[x,y]=RK(0,0.4,2,1)x =0 0.4000 0.8000 1.2000 1.6000 2.0000y =1.0000 0.7328 0.6900 0.8642 1.2248 1.6812（2）2210(1)0y x yxy'⎧=--≤≤⎨-=⎩取步长h=0.1，用四阶标准R-K方法求值。