微分方程常用的两种数值解法：欧拉方法和龙格库塔法

1、 PAGE1 / NUMPAGES33 某师X大学本科毕业论文微分方程常用的两种数值解法：欧拉方法与龙格库塔法学生某XXX院系名称数学与软件科学学院专业名称信息与计算科学班 级2006级 4 班学 号20060640XX指导教师Xxx某师X大学教务处二一年五月微分方程常用的两种数值解法：欧拉方法与龙格库塔法学生某：xxx 指导教师：xx【内容摘要】微分方程是最有生命力的数学分支，在自然科学的许多领域中，都会遇到常微分方程的求解问题。当前计算机的发展为常微分方程的应用及理论研究提供了非常有力的工具，利用计算机解微分方程主要使用数值方法，欧拉方法和龙格库塔方法是求解微分方程最典型常用的数值方法。本

2、文详细研究了这两类数值计算方法的构造过程，分析了它们的优缺点，以及它们的收敛性，相容性，及稳定性。讨论了步长的变化对数值方法的影响和系数不同的同阶龙格库塔方法的差别。通过编制C程序在计算机上实现这两类方法及对一些典型算例的结果分析比较，能更深切体会它们的功能，优缺点及适用场合，从而在实际应用中能对不同类型和不同要求的常微分方程会选取适当的求解方法。关键词：显式单步法 欧拉（Euler）方法 龙格库塔（RungeKutta）方法 截断误差 收敛性Two monly used numerical solution of differential equations：Euler method and

3、 Runge - Kutta methodStudent Name: Xiong Shiying Tutor：Zhang Li【Abstract】The differential equation is the most vitality branch in mathematics. In many domains of natural science, we can meet the ordinary differentialequation solution question. Currently, the development of puter has provided the ext

4、remely powerful tool for the ordinary differential equation application and the fundamental research, the puter solving differential equation mainly uses value method.The Euler method and the RungeKutta method are the most typical monly value method to solve the differential equation. This article d

5、issects the structure process of these two kinds of values monly value method to solve the analyses their good and bad points, to their astringency, the patibility, and the stability has made the proof. At the same time, the article discuss the length of stride to the numerical method changing influ

6、ence and the difference of the coefficient different same step Rungekutta method. Through establishing C program on the puter can realize these two kind of methods, Anglicizing some models of calculate example result can sincerely realize their function, the advantage and disadvantage points and the

7、 suitable situation, thus the suitable solution method can be selected to solve the different type and the different request ordinary differential equation in the practical application .Keywords:Explicit single-step processEuler methodRungeKutta methodtruncation errorconvergence目 录微分方程常用的两种数值解法：欧拉方法

8、与龙格库塔法前言 常微分方程的形成与发展是和力学、天文学、物理学以及其他科学技术的发展密切相关的。数学其他分支的新发展，如复变函数、群、组合拓扑学等，都对常微分方程的发展产生了深刻的影响，当前计算机的发展更是为常微分方程的应用及理论研究提供了非常有力的工具。牛顿在研究天体力学和机械力学的时候，利用了微分方程这个工具，从理论上得到了行星运动规律。后来，法国天文学家勒维烈和英国天文学家亚当斯使用微分方程各自计算出那时尚未发现的海王星的位置。这些都使数学家更加深信微分方程在认识自然、改造自然方面的巨大力量，微分方程也就成了最有生命力的数学分支。然而，我们知道，只有少数十分简单的微分方程能够用初等方法

9、求得它们的解，多数情形只能利用近似方法求解。在常微分方程课中的级数解法，逐步逼近法等就是近似解法。这些方法可以给出解的近似表达式，通常称为近似解析方法。还有一类近似方法称为数值方法，它可以给出解在一些离散点上的近似值，利用计算机解微分方程主要使用数值方法。本文主要讨论一阶常微分方程初值问题 （1.1）在区间上的数值解法，其中为关于，的已知函数，为给定的初始值，将上述问题的精确解记为。该问题常用的数值解法有：欧拉（Euler）方法、龙格库塔（RungeKutta）方法及一些常用的线性多步法。本文重点介绍欧拉（Euler）方法和龙格库塔（RungeKutta）方法。并对这两种方法编制程序，体会它们

10、的功能、优缺点及适用场合，对不同类型常微分方程会选取适当的求解方法。 基本概念和准备知识一阶常微分方程初值问题是：其中是平面上某个区域上的连续函数，式（1.1.1）的微分方程一般有无穷多个解，式（1.1.2）是确定解的初始条件，如果一元函数对一切满足（1）；（2）；（3）存在而且；则称是初值问题（1.1）在区间上的解。误差：假定在计算时，用到的前一步的值是准确的即，把用计算得到的近似值记为，估计误差：= y(xn+1)，这种误差称为局部截断误差。如果不作这一假定，在每一步计算时除局部截断误差以外，还有由于前一步不准确而引起的误差，称为总体截断误差。收敛性：对于解初值问题的数值方法，我们希望它产

11、生的数值解收敛于初值问题的准确解，“收敛性”的一般定义为：对于所有满足引理1.1条件的初值问题（1.1），如果有一种显式单步法：产生的近似解，对于任意固定的，均有，则称该显式单步法是收敛的。相容性：显式单步法（1.2.1）称为与原微分方程相容，如果 （1.2.3）成立。并称式（1.2.3）为相容性条件。稳定性：在实际计算中，一方面初始值不一定精确，往往带有一定误差，同时由于计算机字长有限，在计算过程中有舍入误差，而且这种误差在式（1.2.1）的递推过程中会传递下去，对以后的结果产生影响。因此要考虑舍入误差的积累是否会得到控制，也即要考虑数值方法的稳定性。当时，若舍入误差引起的后果是有限的，则可

12、以认为该方法是数值稳定的。2 欧拉方法2.1欧拉方法简介对常微分方程初值问题（1.1）用数值方法求解时，我们总是认为（1.1）的解存在且唯一。欧拉方法是解初值问题的最简单、最原始的数值方法，它是显式单步法。下面介绍几种导出欧拉法的途径，每个途径皆可以推导出更为有效的数值方法。（1）Taylor展开在点将作Taylor展开得：当充分小时略去误差项，并注意到，得，以近似代替，以近似代替，且用“=”代替“”得差分方程初值问题：， (2.1.1)用式(2.1.1)求解初值问题（1.1）的方法称为欧拉方法。（2）数值微分由导数的定义知，对于充分小的整理得，对此作相似的处理也可以得到欧拉方法(2.1.1)

13、。（3）数值积分在区间对积分得 （2.1.2）用数值积分的左矩形公式计算式（2.1.2）右端的积分，得，于是同样可以得到欧拉方法(2.1.1)。（4）多项式插值利用解和其导数在点的值，作一次埃尔米特插值，得到关于的插值多项式：，用近似代替就得到欧拉方法。（5）待定系数法在第步，已知和，利用这两个值估计出下一步的，将已知的值与估计值作线性组合：，其中， 为待定系数。为确定这两个参数，要求这个估计值对和（为常数）精确成立。如果，则，得到方程：，得。如果，则，这样有：，说明，这样估计值为：，即为欧拉方法。欧拉方法的几何意义：由点斜式得切线方程等步长为，则，可由切线方程算出：，逐步计算出在点的值：，用

14、分段的折线逼近函数为“折线法”而非“切线法”，除第一个点是曲线上的切线，其它都不是。0图1欧拉方法的几何意义2.2欧拉方法的截断误差，收敛性，相容性，稳定性设，把在处展开成泰勒级数，即再由欧拉方法： 两式相减得欧拉方法的局部截断误差为：，若在上充分光滑，且令，则，故欧拉方法是一阶方法或具有一阶精度。欧拉方法的增量函数就是，由引理1.3、引理1.4知当满足Lipschitz条件时欧拉方法是收敛的而且是相容的。用欧拉方法求解典型方程（1.2.4）的计算公式为：，有 。要让，必须有，因此欧拉方法的绝对稳定域为 ：，当为实数时，绝对稳定区间为。在复平面上，是以1为半径、以为圆心的内部。3 龙格库塔法3

15、.1 龙格库塔法的基本思想为了导出龙格库塔法的一般公式，我们取如下的线性组合形式： （3.2.3）其中 （3.2.4）即 ,；a21,a31除c1=0外均为待定系数。显然用公式（3.2.3）每计算一个新值要计算函数的值s次，又因每个都能以一种明显的方式由,计算出来，故将公式（3.2.3）称为s级显式龙格库塔法。s级显式龙格库塔法又可以写成下面既简洁又直观的阵列形式： 03.2 二、三、四级龙格库塔法 常见的二级二阶龙格库塔法有：中点方法（取）二阶Heun方法（取）取s = 3，完全仿照上述方法推导出三阶龙格库塔公式。这时参数满足下列条件这个方程组解也不唯一，从而可以得到不同的三级三阶龙格库塔法

16、。两个较为常见的三级三阶龙格库塔法是：Kutta方法（取=， = 1，a21 = ，a31 = ，a32 = 2， = ， = ， = ）将它代入（3.2.3）得(3.3.2)三阶Heun方法（取=， = ，a21 = ，a31 = 0，a32 = ， = ， = ， =）将它代入（3.2.3）得通常人们所说的龙格库塔法是指四阶而言的。取s=4，我们同样可以仿照二阶的情形推导出此公式。这样就得到如下含12个未知数但仅有10个方程的方程组此方程组的解同样不是唯一的，从而有许多不同的四级四阶龙格库塔法，最常用的两个四阶公式是：经典龙格库塔法 （3.3.3）RKG(Runge-Kutta-Gill)

17、方法 (3.3.4)经典龙格库塔方法的几何解释：将（3.3.3）式中的改记为，从而得到经典龙格库塔方法的计算公式为只考虑区间上的解曲线。由(3.3.3)式可知，曲线在处的斜率，在处的斜率，而和则是在中点处的斜率的两个近似值，如图2：图2另一方面有：，用Simpson积分公式逼近此式右端的积分，得上式右端要计算的三个值。由前面的结果，可取，见图3：图3代入后得到 （3.3.5）另外，从(3.3.3)式容易看出，当与无关时，公式(3.3.5)实际上就是用标准的Simpson积分公式计算积分得到的，因此，可以将经典龙格库塔方法视为是用Simpson积分公式的推广形式计算积分而得到的数值方法。经典龙格

18、库塔法的稳定性：增量函数：，其中代入后得于是有=从而得出龙格库塔方法的绝对稳定区间为 1经典龙格库塔方法的算法框图：开 始读入输出=结 束图4 经典龙格库塔法算法框图4 应用举例建模及其分析4.1 欧拉方法解题及其数学模型4.1.1 问题提出例4.1.1假定某公司的净资产因资产本身产生了利息而以4%的年利率增长，同时，该公司以每年100万的数额支付职工工资。净资产的微分方程为 （t以年为单位）分别以初始值问题：用Euler公式预测公司24后的净资产趋势。4.1.2模型建立 分析：这是求微分方程的数值积分，为的是预测公司24年后的净资产趋势。确定变量：设净资产是时间的微分函数，不妨设变量t为时间

19、（以年为单位）。设为第n年的净资产，为第n+1年的净资产，利息以每年4%的速度增长，且公司每年支付职工工资为100万，则第n+1年的净资产增长数额为，由于已得第n年的净资产为，则第n年的净资产加上第n+1年的净资产增长数额就得到第n+1年的净资产。归纳公式： （4-1）确定其微分方程的标准形式，这就是例4.1.1的微分方程模型。4.1.3 解决问题分别以代入，x表示利息赢利低于工资支出的数额，y表示利息赢利与工资支出平衡的数额，z表示利息赢利高于工资支出的数额，计算结果见表1.表1nn11460.25003540.13834.92625004165.0721418.425003581.6147

20、68.32425004231.6831375.1425003624.8615699.05625004300.9441330.1425003669.8616627.01925004372.9851283.3525003716.6517552.125004447.961234.6825003765.32184748271184.0725003815.93193938581131.4325003868.5720308.87725004691.1291076.6925003923.3121221.23225004778.77101019.762500

21、3980.2422130.08125004869.9211960.54625004039.452335.284525004964.7212898.96825004101.0324-63.304225005063.3从表1可以看到当利息赢利低于工资的支出，公司的净资产逐年减少，以致净资产为负值；当利息赢利与工资的支出平衡时，公司的净资产每年保持不变；当利息赢利超过工资的支出，公司的净资产稳步增长。再用欧拉方法求解，每计算一步，依次需要计算1次、2次和4次函数的值，为了比较在计算量相同的条件下近似解的精度，步长分别取0.025、0.05、0.1。我对欧拉的C程序做了一些优化，加入了计算误差，计算结

22、果如下（表示近似解、表示准确解、表示误差）：This is Euler method:please input a b and y0：0 1 1 please input N：40 x0=0.000000 y0=1.000000 y10=1.000000 e0=0.000000 x1=0.025000 y1=1.025000 y11=1.024695 e1=-0.000305x2=0.050000 y2=1.049405 y12=1.048809 e2=-0.000597x3=0.075000 y3=1.073258 y13=1.072381 e3=-0.000878x4=0.100000 y

23、4=1.096596 y14=1.095445 e4=-0.001151x5=0.125000 y5=1.119451 y15=1.118034 e5=-0.001417x6=0.150000 y6=1.141854 y16=1.140175 e6=-0.001679x7=0.175000 y7=1.163832 y17=1.161895 e7=-0.001937x8=0.200000 y8=1.185410 y18=1.183216 e8=-0.002194x9=0.225000 y9=1.206609 y19=1.204159 e9=-0.002450 x10=0.250000 y10=1

24、.227451 y110=1.224745 e10=-0.002706x11=0.275000 y11=1.247954 y111=1.244990 e11=-0.002964x12=0.300000 y12=1.268134 y112=1.264911 e12=-0.003223x13=0.325000 y13=1.288009 y113=1.284523 e13=-0.003486x14=0.350000 y14=1.307593 y114=1.303841 e14=-0.003753x15=0.375000 y15=1.326900 y115=1.322876 e15=-0.004024

25、x16=0.400000 y16=1.345941 y116=1.341641 e16=-0.004300 x17=0.425000 y17=1.364730 y117=1.360147 e17=-0.004583x18=0.450000 y18=1.383278 y118=1.378405 e18=-0.004873x19=0.475000 y19=1.401594 y119=1.396424 e19=-0.005170 x20=0.500000 y20=1.419689 y120=1.414214 e20=-0.005475x21=0.525000 y21=1.437572 y121=1.

26、431782 e21=-0.005790 x22=0.550000 y22=1.455251 y122=1.449138 e22=-0.006113x23=0.575000 y23=1.472735 y123=1.466288 e23=-0.006447x24=0.600000 y24=1.490032 y124=1.483240 e24=-0.006792x25=0.625000 y25=1.507149 y125=1.500000 e25=-0.007149x26=0.650000 y26=1.524093 y126=1.516575 e26=-0.007518x27=0.675000 y

27、27=1.540872 y127=1.532971 e27=-0.007900 x28=0.700000 y28=1.557490 y128=1.549193 e28=-0.008297x29=0.725000 y29=1.573955 y129=1.565248 e29=-0.008708x30=0.750000 y30=1.590273 y130=1.581139 e30=-0.009134x31=0.775000 y31=1.606449 y131=1.596872 e31=-0.009577x32=0.800000 y32=1.622489 y132=1.612452 e32=-0.0

28、10037x33=0.825000 y33=1.638397 y133=1.627882 e33=-0.010515x34=0.850000 y34=1.654180 y134=1.643168 e34=-0.011013x35=0.875000 y35=1.669842 y135=1.658312 e35=-0.011530 x36=0.900000 y36=1.685388 y136=1.673320 e36=-0.012068x37=0.925000 y37=1.700823 y137=1.688194 e37=-0.012629x38=0.950000 y38=1.716151 y13

29、8=1.702939 e38=-0.013212x39=0.975000 y39=1.731376 y139=1.717556 e39=-0.013820 x40=1.000000 y40=1.746504 y140=1.732051 e40=-0.014453this is improve Euler method please input a b and y0:：0 1 1 please input N：20 x0=0.000000 y0=1.000000 y10=1.000000 e0=0.000000 x1=0.050000 y1=1.048869 y11=1.048809 e1=-0

30、.000060 x2=0.100000 y2=1.095561 y12=1.095445 e2=-0.000116x3=0.150000 y3=1.140345 y13=1.140175 e3=-0.000169x4=0.200000 y4=1.183437 y14=1.183216 e4=-0.000221x5=0.250000 y5=1.225017 y15=1.224745 e5=-0.000272x6=0.300000 y6=1.265236 y16=1.264911 e6=-0.000325x7=0.350000 y7=1.304219 y17=1.303841 e7=-0.0003

31、78x8=0.400000 y8=1.342074 y18=1.341641 e8=-0.000433x9=0.450000 y9=1.378896 y19=1.378405 e9=-0.000492x10=0.500000 y10=1.414766 y110=1.414214 e10=-0.000553x11=0.550000 y11=1.449756 y111=1.449138 e11=-0.000618x12=0.600000 y12=1.483927 y112=1.483240 e12=-0.000687x13=0.650000 y13=1.517337 y113=1.516575 e

32、13=-0.000762x14=0.700000 y14=1.550035 y114=1.549193 e14=-0.000841x15=0.750000 y15=1.582066 y115=1.581139 e15=-0.000927x16=0.800000 y16=1.613472 y116=1.612452 e16=-0.001020 x17=0.850000 y17=1.644289 y117=1.643168 e17=-0.001121x18=0.900000 y18=1.674551 y118=1.673320 e18=-0.001230 x19=0.950000 y19=1.70

33、4288 y119=1.702939 e19=-0.001349x20=1.000000 y20=1.733529 y120=1.732051 e20=-0.0014784.1.4 结论这道题用普通的微分方程也能列式求解，关键在于如何预测若干年后的净资产趋势，用普通的微分方程就无法进行预测，且人工计算量相当大，这里我使用欧拉方法可以计算出精度以及误差，通过电脑运行程序就可以预测出若干年后的净资产情况，欧拉方法的使用使得解题更加方便且精确。4.2 龙格库塔解题及其数学模型4.2.1 问题提出 例4.2.1 两种果树寄生虫，其数量分别是其中一种寄生虫以吃另一种寄生虫为生，两种寄生虫的增长函数如下列

34、常微分方程组所示：问题：预测3年后这一对寄生虫的数量。4.2.2 模型建立分析：这是一个典型的常微分方程例题，要求预测出3年后这对寄生虫的数量。确定变量：假定时间是寄生虫数量的积分函数，不妨设变量t为时间（以年为单位）。由题知有两种寄生虫u和v，u寄生虫年增长函数为，v寄生虫年增长函数为，初始值：u寄生虫为1.6，v寄生虫为1.2，由于其中一种寄生虫以吃另一种寄生虫为生，我们可建立u和v的关联函数f(u,v),g(u,v)。归纳后得到公式： （4-2）（4-2）即为例4.2.1所述问题的微分方程模型。4.2.3 解决问题在本题中用Euler预估校正公式取h=1，计算结果如表2.表2t/年u(t

35、)v(t)11.61.221.024571.2683430.6409121.336640.3912111.41077我把龙格库塔的C程序进行编辑后得到结果:（ei表示误差）please input a b and y0:：0 1 1 please input N： 10 x0=0.000000 y0=1.000000 y10=1.000000 e0=0.000000 x1=0.100000 y1=1.095446 y11=1.095445 e1=-0.000000 x2=0.200000 y2=1.183217 y12=1.183216 e2=-0.000001x3=0.300000 y3=1

36、.264912 y13=1.264911 e3=-0.000001x4=0.400000 y4=1.341642 y14=1.341641 e4=-0.000001x5=0.500000 y5=1.414215 y15=1.414214 e5=-0.000002x6=0.600000 y6=1.483242 y16=1.483240 e6=-0.000002x7=0.700000 y7=1.549196 y17=1.549193 e7=-0.000003x8=0.800000 y8=1.612455 y18=1.612452 e8=-0.000004x9=0.900000 y9=1.67332

37、4 y19=1.673320 e9=-0.000004x10=1.000000 y10=1.732056 y110=1.732051 e10=-0.0000054.2.3 结论从上面的结果可以分析，用每一种方法计算节点=0.1、0.20.9、1.0，上的值都需要计算4次的值，即它们的计算量基本相同，其结果是经典的龙格库塔方法的精度最好，龙格库塔方法的推导是基于Taylor级数方法的，因而在使用高阶龙格库塔方法计算时可以很精确的推算出寄生虫每一年的数量。4.3 分别使用二阶、三阶龙格库塔法解初值问题对一些特殊的微分方程，使用欧拉方法和低阶的龙格库塔方法也能达到很高的精度，例如：微分方程的解析解是

38、一次函数，则用欧拉方法求得的数值解与准确解相符，微分方程的解析解是二次函数，则用二阶龙格库塔方法求得的数值解与准确解相符。微分方程的解析解是三次多项式，用三阶龙格库塔方法求得的数值解与准确解相符。4.3.1 建立模型建立初值问题模型： (4.1)4.3.2 用不同的龙格库塔法解（4.1）初值问题模型 （1）用二阶龙格库塔方法求解初值问题（4.1）结果如下：this is the second-order Runge-Kutta(Heun) method please input a b and y0： 0 1 1 please input N：10 x0=0.000000 y0=1.00000

39、0 y10=1.000000 e0=0.000000 x1=0.100000 y1=1.110000 y11=1.110000 e1=0.000000 x2=0.200000 y2=1.240000 y12=1.240000 e2=0.000000 x3=0.300000 y3=1.390000 y13=1.390000 e3=0.000000 x4=0.400000 y4=1.560000 y14=1.560000 e4=0.000000 x5=0.500000 y5=1.750000 y15=1.750000 e5=0.000000 x6=0.600000 y6=1.960000 y16=

40、1.960000 e6=0.000000 x7=0.700000 y7=2.190000 y17=2.190000 e7=0.000000 x8=0.800000 y8=2.440000 y18=2.440000 e8=0.000000 x9=0.900000 y9=2.710000 y19=2.710000 e9=0.000000 x10=1.000000 y10=3.000000 y110=3.000000 e10=0.000000this is improve Euler methodplease input a b and y0：0 1 1 please input N：10 x0=0

41、.000000 y0=1.000000 y10=1.000000 e0=0.000000 x1=0.100000 y1=1.110000 y11=1.110000 e1=0.000000 x2=0.200000 y2=1.240000 y12=1.240000 e2=0.000000 x3=0.300000 y3=1.390000 y13=1.390000 e3=0.000000 x4=0.400000 y4=1.560000 y14=1.560000 e4=0.000000 x5=0.500000 y5=1.750000 y15=1.750000 e5=0.000000 x6=0.60000

42、0 y6=1.960000 y16=1.960000 e6=0.000000 x7=0.700000 y7=2.190000 y17=2.190000 e7=0.000000 x8=0.800000 y8=2.440000 y18=2.440000 e8=0.000000 x9=0.900000 y9=2.710000 y19=2.710000 e9=0.000000 x10=1.000000 y10=3.000000 y110=3.000000 e10=0.000000this is the second-order Runge-Kutta(middle) methodplease inpu

43、t a b and y0 ：0 1 1 please input N：10 x0=0.000000 y0=1.000000 y10=1.000000 e0=0.000000 x1=0.100000 y1=1.110000 y11=1.110000 e1=0.000000 x2=0.200000 y2=1.240000 y12=1.240000 e2=0.000000 x3=0.300000 y3=1.390000 y13=1.390000 e3=0.000000 x4=0.400000 y4=1.560000 y14=1.560000 e4=0.000000 x5=0.500000 y5=1.

44、750000 y15=1.750000 e5=0.000000 x6=0.600000 y6=1.960000 y16=1.960000 e6=0.000000 x7=0.700000 y7=2.190000 y17=2.190000 e7=0.000000 x8=0.800000 y8=2.440000 y18=2.440000 e8=0.000000 x9=0.900000 y9=2.710000 y19=2.710000 e9=0.000000 x10=1.000000 y10=3.000000 y110=3.000000 e10=0.000000（2）用三阶龙格库塔方法求解初值问题模型

45、： （4.2）结果如下：this is the third-order Heun method:please input a b and y0：0 2 1 please input N：10 x0=0.000000 y0=1.000000 y10=1.000000 e0=0.000000 x1=0.200000 y1=1.248000 y11=1.248000 e1=0.000000 x2=0.400000 y2=1.624000 y12=1.624000 e2=0.000000 x3=0.600000 y3=2.176000 y13=2.176000 e3=0.000000 x4=0.800

46、000 y4=2.952000 y14=2.952000 e4=0.000000 x5=1.000000 y5=4.000000 y15=4.000000 e5=-0.000000 x6=1.200000 y6=5.368000 y16=5.368000 e6=0.000000 x7=1.400000 y7=7.104000 y17=7.104000 e7=-0.000000 x8=1.600000 y8=9.256000 y18=9.256000 e8=-0.000000 x9=1.800000 y9=11.872000 y19=11.872000 e9=-0.000000 x10=2.00

47、0000 y10=15.000000 y110=15.000000 e10=-0.000000this is kutta methodplease input a b and y0：0 2 1 please input N：10 x0=0.000000 y0=1.000000 y10=1.000000 e0=0.000000 x1=0.200000 y1=1.248000 y11=1.248000 e1=0.000000 x2=0.400000 y2=1.624000 y12=1.624000 e2=0.000000 x3=0.600000 y3=2.176000 y13=2.176000 e

48、3=0.000000 x4=0.800000 y4=2.952000 y14=2.952000 e4=0.000000 x5=1.000000 y5=4.000000 y15=4.000000 e5=0.000000 x6=1.200000 y6=5.368000 y16=5.368001 e6=0.000000 x7=1.400000 y7=7.104000 y17=7.104000 e7=-0.000000 x8=1.600000 y8=9.256001 y18=9.256001 e8=0.000000 x9=1.800000 y9=11.872001 y19=11.872001 e9=0

49、.000000 x10=2.000000 y10=15.000001 y110=15.000000 e10=-0.0000014.3.3 结论如果常微分方程的解析解是任意的次多项式，总存在阶龙格库塔方法，用它求解该微分方程的所得的数值解与准确解相符。证明：龙格库塔方法的推导是基于Taylor级数方法的，阶龙格库塔方法的局部截断误差为：，其中：，为常数，如果常微分方程的解析解是任意的次多项式，则有，因此局部截断误差等于零。从而上述结论成立。5 步长算法模型及其分析下面我们分析步长对不同数值方法的影响。用欧拉方法解初值问题建立初值模型： （4.3）当h=0.2时this is Euler meth

50、odplease input a b and y0：0 1 1 please input N： 5x0=0.000000 y0=1.000000 y10=1.000000 e0=0.000000 x1=0.200000 y1=1.200000 y11=1.183216 e1=-0.016784x2=0.400000 y2=1.373333 y12=1.341641 e2=-0.031693x3=0.600000 y3=1.531495 y13=1.483240 e3=-0.048256x4=0.800000 y4=1.681085 y14=1.612452 e4=-0.068633x5=1.0

51、00000 y5=1.826948 y15=1.732051 e5=-0.094898当h=0.1时this is Euler method :please input a b and y0：0 1 1 please input N： 10 x0=0.000000 y0=1.000000 y10=1.000000 e0=0.000000 x1=0.100000 y1=1.100000 y11=1.095445 e1=-0.004555x2=0.200000 y2=1.191818 y12=1.183216 e2=-0.008602x3=0.300000 y3=1.277438 y13=1.26

52、4911 e3=-0.012527x4=0.400000 y4=1.358213 y14=1.341641 e4=-0.016572x5=0.500000 y5=1.435133 y15=1.414214 e5=-0.020920 x6=0.600000 y6=1.508966 y16=1.483240 e6=-0.025727x7=0.700000 y7=1.580338 y17=1.549193 e7=-0.031145x8=0.800000 y8=1.649784 y18=1.612452 e8=-0.037332x9=0.900000 y9=1.717780 y19=1.673320

53、e9=-0.044460 x10=1.000000 y10=1.784771 y110=1.732051 e10=-0.052720用改进的欧拉方法解初值问题模型（4.3）。当h=0.2时this is improve Euler method :please input a b and y0：0 1 1 please input N：5x0=0.000000 y0=1.000000 y10=1.000000 e0=0.000000 x1=0.200000 y1=1.186667 y11=1.183216 e1=-0.003451x2=0.400000 y2=1.348312 y12=1.34

54、1641 e2=-0.006672x3=0.600000 y3=1.493704 y13=1.483240 e3=-0.010464x4=0.800000 y4=1.627861 y14=1.612452 e4=-0.015410 x5=1.000000 y5=1.754205 y15=1.732051 e5=-0.022154当h=0.1时this is improve Euler method :please input a b and y0：0 1 1 please input N： 10 x0=0.000000 y0=1.000000 y10=1.000000 e0=0.000000

55、x1=0.100000 y1=1.095909 y11=1.095445 e1=-0.000464x2=0.200000 y2=1.184097 y12=1.183216 e2=-0.000881x3=0.300000 y3=1.266201 y13=1.264911 e3=-0.001290 x4=0.400000 y4=1.343360 y14=1.341641 e4=-0.001719x5=0.500000 y5=1.416402 y15=1.414214 e5=-0.002188x6=0.600000 y6=1.485955 y16=1.483240 e6=-0.002715x7=0.

56、700000 y7=1.552514 y17=1.549193 e7=-0.003320 x8=0.800000 y8=1.616474 y18=1.612452 e8=-0.004023x9=0.900000 y9=1.678166 y19=1.673320 e9=-0.004846x10=1.000000 y10=1.737867 y110=1.732051 e10=-0.005816用二级二阶龙格库塔方法（二阶Heun方法）解初值问题模型（4.3）。当h=0.2时this is the second-order Runge-Kutta method :please input a b a

57、nd y0：0 1 1 please input N：5x0=0.000000 y0=1.000000 y10=1.000000 e0=0.000000 x1=0.200000 y1=1.184706 y11=1.183216 e1=-0.001490 x2=0.400000 y2=1.344644 y12=1.341641 e2=-0.003003x3=0.600000 y3=1.488062 y13=1.483240 e3=-0.004822x4=0.800000 y4=1.619653 y14=1.612452 e4=-0.007201x5=1.000000 y5=1.742496 y1

58、5=1.732051 e5=-0.010445当h=0.1时this is the second-order Runge-Kutta method :please input a b and y0：0 1 1 please input N：10 x0=0.000000 y0=1.000000 y10=1.000000 e0=0.000000 x1=0.100000 y1=1.095625 y11=1.095445 e1=-0.000180 x2=0.200000 y2=1.183572 y12=1.183216 e2=-0.000356x3=0.300000 y3=1.265449 y13=1

59、.264911 e3=-0.000538x4=0.400000 y4=1.342374 y14=1.341641 e4=-0.000733x5=0.500000 y5=1.415162 y15=1.414214 e5=-0.000948x6=0.600000 y6=1.484431 y16=1.483240 e6=-0.001191x7=0.700000 y7=1.550663 y17=1.549193 e7=-0.001470 x8=0.800000 y8=1.614246 y18=1.612452 e8=-0.001794x9=0.900000 y9=1.675493 y19=1.6733

60、20 e9=-0.002173x10=1.000000 y10=1.734671 y110=1.732051 e10=-0.002620用三级三阶龙格库塔方法（Kutta方法）解初值问题模型（4.3）。当h=0.2时this is kutta method :please input a b and y0 ：0 1 1 please input N：5x0=0.000000 y0=1.000000 y10=1.000000 e0=0.000000 x1=0.200000 y1=1.183244 y11=1.183216 e1=-0.000028x2=0.400000 y2=1.341729 y