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

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

2、构造过程，分析了它们的优缺点，以及它们的收敛性，相容性，及稳定性。讨论了步长的变化对数值方法的影响和系数不同的同阶龙格一库塔方法的差别。通过编制C程序在计算机上实现这两类方法及对一些典型算例的结果分析比较，能更深切体会它们的功能，优缺点及适用场合，从而在实际应用中能对不同类型和不同要求的常微分方程会选取适当的求解方法。关键词：显式单步法欧拉(Euler)方法龙格一库塔(Runge-Kutta)方法截断误差收敛性Twocommonlyusednumericalsolutionofdifferentialequations:EulermethodandRunge-Kuttamethod:Zhang

3、LiStudentName:XiongShiyingTutorAbstractThedifferentialequationisthemostvitalitybranchinmathematics.Inmanydomainsofnaturalscience,wecanmeettheordinarydifferentialequationsolutionquestion.Currently,thedevelopmentofcomputerhasprovidedtheextremelypowerfultoolfortheordinarydifferentialequationapplication

4、andthefundamentalresearch,thecomputersolvingdifferentialequationmainlyusesvaluemethod.TheEulermethodandtheRungeKuttamethodarethemosttypicalcommonlyvaluemethodtosolvethedifferentialequation.Thisarticledissectsthestructureprocessofthesetwokindsofvaluescommonlyvaluemethodtosolvetheanalysestheirgoodandb

5、adpoints,totheirastringency,thecompatibility,andthestabilityhasmadetheproof.Atthesametime,thearticlediscussthelengthofstridetothenumericalmethodchanginginfluenceandthedifferenceofthecoefficientdifferentsamestepRungekuttamethod.ThroughestablishingCprogramonthecomputercanrealizethesetwokindofmethods,A

6、nglicizingsomemodelsofcalculateexampleresultcansincerelyrealizetheirfunction,theadvantageanddisadvantagepointsandthesuitablesituation,thusthesuitablesolutionmethodcanbeselectedtosolvethedifferenttypeandthedifferentrequestordinarydifferentialequationinthepracticalapplication.Keywords：Explicitsingle-s

7、tepprocessEulermethodRungeKuttamethodtruncationerrorconvergence微分方程常用的两种数值解法：欧拉方法与龙格一库塔法、/、一刖百常微分方程的形成与发展是和力学、天文学、物理学以及其他科学技术的发展密切相关的。数学其他分支的新发展，如复变函数、群、组合拓扑学等，都对常微分方程的发展产生了深刻的影响，当前计算机的发展更是为常微分方程的应用及理论研究提供了非常有力的工具。牛顿在研究大体力学和机械力学的时候，利用了微分方程这个工具，从理论上得到了行星运动规律。后来，法国天文学家勒维烈和英国天文学家亚当斯使用微分方程各自计算出那时尚未发现的海王

8、星的位置。这些都使数学家更加深信微分方程在认识自然、改造自然方面的巨大力量，微分方程也就成了最有生命力的数学分支。然而，我们知道，只有少数十分简单的微分方程能够用初等方法求得它们的解，多数情形只能利用近似方法求解。在常微分方程课中的级数解法，逐步逼近法等就是近似解法。这些方法可以给出解的近似表达式，通常称为近似解析方法。还有一类近似方法称为数值方法，它可以给出解在一些离散点上的近似值，利用计算机解微分方程主要使用数值方法。本文主要讨论一阶常微分方程初值问题(1.1)f(x,y)y(x0)=y0在区间a,b上的数值解法，其中f(x,y)为关于x,y的已知函数，y。为给定的初始值，将上述问题的精确

9、解记为y(x)0该问题常用的数值解法有：欧拉(Euler)方法、龙格一库塔(Runge-Kutta)方法及一些常用的线性多步法。本文重点介绍欧拉(Euler)方法和龙格一库塔(Runge-Kutta)方法。并对这两种方法编制程序，体会它们的功能、优缺点及适用场合，对不同类型常微分方程会选取适当的求解方法。1基本概念和准备知识一阶常微分方程初值问题是:(1.1.1)(1.1.2)dy=f(x,y)Sdxj(x。)=y。其中f(x,y)是平面上某个区域D上的连续函数，式(1.1.1)的微分方程一般有无穷多个解，式(1.1.2)是确定解的初始条件，如果一元函数y(x)对一切aExMb满足(1) (x

10、,y(x”D;(2) y(x0)=y0;(3) y存在而且y(x)=f(x,y(x)；则称y(x)是初值问题(1.1)在区间a,b上的解。误差：假定在计算时，用到的前一步的值yn是准确的(即yn=y(xn),把用y(xn)计算得到的近似值记为14估计误差：En+=y(xn+l)-1书，这种误差称为局部截断误差。如果不作这一假定，在每一步计算时除局部截断误差以外，还有由于前一步不准确而引起的误差，称为总体截断误差。收敛性：对于解初值问题的数值方法，我们希望它产生的数值解收敛于初值问题的准确解，“收敛性”的一般定义为：对于所有满足引理1.1条件的初值问题(1.1),如果有一种显式单步法：Yn1=Y

11、nh；：(xn,yn,h)产生的近似解，对于任意固定的xn=x0+nh,均有limyn=y(xn),则称该显式单步h_0法是收敛的。相容性：显式单步法(1.2.1)称为与原微分方程相容，如果巴x,y(x),0)=f(x,y(x)(1.2.3)成立。并称式(1.2.3)为相容性条件。稳定性：在实际计算中，一方面初始值y。不一定精确，往往带有一定误差，同时由于计算机字长有限，在计算过程中有舍入误差，而且这种误差在式(1.2.1)的递推过程中会传递下去，对以后的结果产生影响。因此要考虑舍入误差的积累是否会得到控制，也即要考虑数值方法的稳定性。当nT8时，若舍入误差引起的后果是有限的，则可以认为该方法

12、是数值稳定的。2欧拉方法2.1欧拉方法简介对常微分方程初值问题(1.1)用数值方法求解时，我们总是认为(1.1)的解存在且唯一。欧拉方法是解初值问题的最简单、最原始的数值方法，它是显式单步法。下面介绍几种导出欧拉法的途径，每个途径皆可以推导出更为有效的数值方法。(1)Taylor展开在x-点将y(xn+)=y(xn+h)作Taylor展开得:y(Xn1)=y(Xnh)h2-n(xn，Xn.1)=y(Xn)hy(Xn)y(；n)h2当h充分小时略去块差项一y(sn)并注息到y(Xn)=f(Xn,y(Xn),得2y(XnQ定y(Xn)+hf(Xn,y(Xn),以丫口近似代替y(Xn),以丫口书近似

13、代替y(Xn中)，且用“二”代替得差分方程初值问题：In书=Vn+hf(Xn,yn)ba,_.,yo=y(Xo)Xn=X0+nh,n=0,1,2,，N1,N=b-(2.1.1)、h用式(2.1.1)求解初值问题(1.1)的方法称为欧拉方法。(2)数值微分由导数的定义知，对于充分小的hy(Xn1)-y(Xn)y(Xn)=f(Xn,y(Xn)h整理得y(Xn书)定y(Xn)+hf(Xn,y(Xn),对此作相似的处理也可以得到欧拉方法(2.1.(1)(3)数值积分在区间Xn，Xn书对y=f(X,y(X)积分得Xn1y(Xn+)=y(Xn)十f(X,y(X)dX(2.1.2)Xn用数值积分的左矩形公式

14、计算式(2.1.2)右端的积分，得y(Xn+)定y(Xn)十hf(Xn,y(Xn),于是同样可以得到欧拉方法(2.1.1)。(4)多项式插值利用解和其导数f(X,y)在Xn点的值yn,f(Xn,yn)作一次埃尔米特插值，得到关于h的插值多项式：p(h)=yn+hf(Xn，yn),用p(h)近似代替yn书就得到欧拉方法。(5)待定系数法在第n步，已知yn和yn=f(Xn,yn),利用这两个值估计出下一步的yn书，将已知的值与估计值作线性组合：yn=ayn+Pyn,其中a,P为待定系数。为确定这两个参数，要求这个估计值对y(X)三c和y(X)=cx(c为常数)精确成立。如果这样有:y(x)三c,y

15、(x)=0,得到方程：c=uc+P0,得以=1。如果y(x)=cx，贝Uy(x)=c,CXn+=Otcxn+Pc=c(xn+P),说明P=xn+-xn=h,这样估计值为:yn+=yn+hyn=yn+hf(xn,yn)，即为欧拉方法。欧拉方法的几何意义:由点斜式得切线方程xc,y0.门y。(xx0)f(x。，y。)dyy=y0(x-x)dx等步长为h,则xi-x0=h,可由切线方程算出yi:yi=y0+hf(x0,y0),逐步计算出y=y(x)在xn书点的值：yn=yn+hf(xn,yn),n=0,1,2,，用分段的折线逼近函数为“折线法”而非“切线法”，除第一个点是曲线上的切线，其它都不是。图

16、1欧拉方法的几何意义2.2欧拉方法的截断误差，收敛性，相容性，稳定性设yn=y(xn),把y(xn书)在xn处展开成泰勒级数，即h2.y(xni)y(xn)hy(xn)y()；n*n,Xni2!再由欧拉方法：yni=ynhf(xn,yn)=y(xn)hf(xn,y(xn)=y(xn)hy(xn)两式相减得欧拉方法的局部截断误差为：En书=h2y(君n)，若y(x)在a,b上充分2，一一.1cc光滑，且令M=maxy(x),则En中|-hM=h),故欧拉方法是一阶方法或具有一阶精度。欧拉方法的增量函数中(x,y,h)就是f,由引理1.3、引理1.4知当f满足Lipschitz条件时欧拉方法是收敛

17、的而且是相容的。用欧拉方法求解典型方程(1.2.4)的计算公式为：yn邛=yn十九hyn=(1十九h)yn,有15n斗=(1+?,_h)6n。要让回1同，必须有1+Kh1,因此欧拉方法的绝对稳定域为：1+九h1,当九为实数时，绝对稳定区间为(-2,0)。在复平面上，1十九h1是以1为半径、以-1为圆心的内部。3龙格一库塔法3.1龙格一库塔法的基本思想为了导出龙格一库塔法的一般公式，我们取如下的线性组合形式：syn+=yn+hbiki(3.2.3)i1i1其中ki=f(xn+cih,yn+hZajkj)(3.2.4)jwk=f%,yn)口k2=f(xn+C2h,yn+ha21k1)k3=f(xn

18、C3h,ynh%1kha32k2)1 .b1，b2,，bs；C1,C2,C3,，Cs;a21,a31,，ass除c=0外均为待定系数。显然用公式(3.2.3)每计算一个新值yn书要计算函数f的值s次，又因每个kj都能以一种明显的方式由k1,k2,，kj计算出来，故将公式(3.2.3)称为s级显式龙格一库塔法。s级显式龙格一库塔法又可以写成下面既简洁又直观的阵列形式：c2a21c3a31a32aamcsas1as2assb1b2bsbs3.2二、三、四级龙格一库塔法常见的二级二阶龙格一库塔法有：1中点方法(取b2=1,b1=0,a21=c2=)2n由=yn+hk2，ki=f(Xn,yn)hhk2

19、=f(Xn+-,yn+-ki)312阶Heun方法(取b2=,b1=,a21=c2=)443h-yn噌=yn十(k1+3k4)4,k1=f(Xn,yn)22k2=f(Xn+-h,yn+-hk1)l.33取s=3,完全仿照上述方法推导出三阶龙格一库塔公式。这时参数满足下列条件b+b2+th=1c2b2+c3b3=12I,2-21b2c2+b3c3=3-1b3c2a32-6这个方程组解也不唯一，从而可以得到不同的三级三阶龙格一库塔法。两个较为常见的三级三阶龙格一库塔法是：1Kutta方法(取cz=,c3=1,a21=a3141、,b3=)66将它代入(3.2.3)得h-yn1二yn(k14k2-k

20、3)6k1=f(Xn，yn)k2=f(xnk3=f(Xnhk1，ynh力22h,yn-hk12hk2)(3.3.2)三阶Heun方法（取一，a21=二，322a31-0,a32,3，1，b1一，b2430,b3=3)将它代入(3.2.3)4yn1=ynh(-k1k3)44k1=f(Xn,yn)k2=f(Xnk3=f(XnT，yn冷3322-h,yn-hk2)33通常人们所说的龙格一库塔法是指四阶而言的。取s=4,我们同样可以仿照二阶的情形推导出此公式。这样就得到如下含12个未知数但仅有10个方程的方程组,2,2,2b2c2b3c3b4c42131b3c2a32b4(c2a42c3a43).2.

21、,2b3c2a32b4(c2a42c2a43)=12a31+a32c3a41+a42+a43=c4b1+b2+b3+b4=1，1b2c2b3c2b4c4=一,3,.3,.3b2c2+b35+b4c41b3c2c3a32b4c4(c2a42c3a43)=二81b4c2a32a43二一此方程组的解同样不是唯一的，从而有许多不同的四级四阶龙格一库塔法，最常用的两个四阶公式是：经典龙格一库塔法h.一y5=yn+7(kl+2k2+2k3+k4)6ki=f(Xn,yn)1 ,1.h.、4 k2=f(Xn+h,yn+ki)1hk3=f(Xn十一h,yn+k2)22k4二f(Xnh,ynhk3)(3.3.3)

22、RKG(Runge-Kutta-Gill)方法hi-(-yn+=yn+hk1+(2-2)k2+(2+“；2冰3+k46K=f(Xn,yn)一,h4k2=f(Xn+-,ynhk3=f(Xn+-,ynk4=f(Xn+h,yn2-12-2hk2)hk2hk3)(3.3.4)经典龙格一库塔方法的几何解释:将(3.3.3)式中的kj改记为fj,从而得到经典龙格一库塔方法的计算公式为hyn1=yn,9(fl,f2f3f4)只考虑区间X0,X上的解曲线y=y(X)o由(3.3.3)式可知，曲线y=y(x)在X。处的斜率y(Xo)=f(X。，y(X。)=f1，在为=x。+h处的斜率y(X1)=f(x1,y(X

23、1)刘f4，而f2和f3则是在中点X1：2=X。+h处的斜率y(X1，2)=f(X12,y(X1：2)的两个近似值，如2图2:图2Xi另一方面有：y(x1)-y(x0)=ff(x,y(x)dx,用Simpson积分公式谩近此式右Xo端的积分，得h_一一y(xi)-y(xo)f(xo,y(xo)4f(xi2,y(xi2)f(xi,y(xi)6上式右端要计算f(x,y)的三个值。由前面的结果，可取f(x0,y(x0)=f1，if(xiy(xi)f4,勺fjfj,见图3:代入后得到h4(f2f3)h,、yi=Yo+fl+f4=y0+(fl+2f2+2f3+f4)(3.3.5)626另外，从(3.3.

24、3)式容易看出，当f(x,y)与y无关时，公式(3.3.5)实际上就是用标Xi准的Simpson积分公式计算积分ff(x,y)dx得到的，因此，可以将经典龙格库x0塔方法视为是用Simpson积分公式的推广形式计算积分广+f(x,y)dx而得到的数值Xn方法。经典龙格一库塔法的稳定性：1.增重函数：中(x,y,h)=-(k+2k2+2k3+k4),其中6k1二1yI12k2=1y一1hy2,12,13,2k3-yhyhy24k4-y+；2hy1-3h2y4h3y224代入后得平(x,y,h)=九y+九2hy+1九3h2y+-九3h2y+九4h3y2624于是有-_，21，3h21，3h21，4

25、h3-.hhhh:y2624从而得出龙格一库塔方法的绝对稳定区间为九十九2h十儿3h2十儿3h2十,九4h3|12624-2.785h0经典龙格一库塔方法的算法框图：图4经典龙格一库塔法算法框图4应用举例建模及其分析4.1 欧拉方法解题及其数学模型4.1.1 问题提出例4.1.1假定某公司的净资产因资产本身产生了利息而以4%勺年利率增长，同时，该公司以每年100万的数额支付职工工资。净资产的微分方程为dw=0.04w100（t以年为单位）dt分别以初始值x(0)=1500,y(0)=2500,z(0)=3500,问题：用Euler公式预测公司24后的净资产趋势4.1.2 模型建立分析：这是求微

26、分方程的数值积分，为的是预测公司24年后的净资产趋势。确定变量：设净资产是时间的微分函数，不妨设变量t为时间(以年为单位)。设wn为第n年的净资产，wn+1为第n+1年的净资产，利息以每年4%勺速度增长，且公司每年支付职工工资为100万，则第n+1年的净资产增长数额为(0.04wn100),由于已得第n年的净资产为wn,则第n年的净资产加上第n+1年的净资产增长数额就得到第n+1年的净资产wn+10归纳公式：wn44=wn+h(0.04wn100)=1.04wn100,h=1.(4-1)确定其微分方程的标准形式，这就是例4.1.1的微分方程模型。4.1.3 解决问题w。分别以x0=1500,y

27、。=2500,z0=3500代入，x表示利息赢利低于工资支出的数额，y表示利息赢利与工资支出平衡的数额，z表示利息赢利高于工资支出的数额，计算结果见表1.表1nxnynznnxnynzn11460.25003540.13834.92625004165.0721418.425003581.614768.32425004231.6831375.1425003624.8615699.05625004300.9441330.1425003669.8616627.01925004372.9851283.3525003716.6517552.125004447.961234.6825003765.3218

28、4748271184.0725003815.93193938581131.4325003868.5720308.87725004691.1291076.6925003923.3121221.23225004778.77101019.7625003980.2422130.08125004869.9211960.54625004039.452335.284525004964.7212898.96825004101.0324-63.304225005063.3从表1可以看到当利息赢利低于工资的支出，公司的净资产逐年减少，以致净资产为负值；当利息赢利

29、与工资的支出平衡时，公司的净资产每年保持不变；当利息赢利超过工资的支出，公司的净资产稳步增长。再用欧拉方法求解，每计算一步y,依次需要计算1次、2次和4次函数f的值，为了比较在计算量相同的条件下近似解的精度，步长分别取0.025、0.05、0.1。我对欧拉的C程序做了一些优化，加入了计算误差，计算结果如下（yi表示近似解、y1i表示准确解、ei表示误差）：ThisisEulermethod:pleaseinputabandy0:011pleaseinputN:40x0=0.000000y0=1.000000y10=1.000000e0=0.000000x1=0.025000y1=1.02500

30、0y11=1.024695e1=-0.000305x2=0.050000y2=1.049405y12=1.048809e2=-0.000597x3=0.075000y3=1.073258y13=1.072381e3=-0.000878x4=0.100000y4=1.096596y14=1.095445e4=-0.001151x5=0.125000y5=1.119451y15=1.118034e5=-0.001417x6=0.150000y6=1.141854y16=1.140175e6=-0.001679x7=0.175000y7=1.163832y17=1.161895e7=-0.00193

31、7x8=0.200000y8=1.185410y18=1.183216e8=-0.002194x9=0.225000y9=1.206609y19=1.204159e9=-0.002450x10=0.250000y10=1.227451y110=1.224745e10=-0.002706x11=0.275000y11=1.247954y111=1.244990e11=-0.002964x12=0.300000y12=1.268134y112=1.264911e12=-0.003223x13=0.325000y13=1.288009y113=1.284523e13=-0.003486x14=0.3

32、50000y14=1.307593y114=1.303841e14=-0.003753x15=0.375000y15=1.326900y115=1.322876e15=-0.004024x16=0.400000y16=1.345941y116=1.341641e16=-0.004300x17=0.425000y17=1.364730y117=1.360147e17=-0.004583x18=0.450000y18=1.383278y118=1.378405e18=-0.004873x19=0.475000y19=1.401594y119=1.396424e19=-0.005170x20=0.5

33、00000y20=1.419689y120=1.414214e20=-0.005475x21=0.525000y21=1.437572y121=1.431782e21=-0.005790x22=0.550000y22=1.455251y122=1.449138e22=-0.006113x23=0.575000y23=1.472735y123=1.466288e23=-0.006447x24=0.600000y24=1.490032y124=1.483240e24=-0.006792x25=0.625000y25=1.507149y125=1.500000e25=-0.007149x26=0.6

34、50000y26=1.524093y126=1.516575e26=-0.007518x27=0.675000y27=1.540872y127=1.532971e27=-0.007900x28=0.700000y28=1.557490y128=1.549193e28=-0.008297x29=0.725000y29=1.573955y129=1.565248e29=-0.008708x30=0.750000y30=1.590273y130=1.581139e30=-0.009134x31=0.775000y31=1.606449y131=1.596872e31=-0.009577x32=0.8

35、00000y32=1.622489y132=1.612452e32=-0.010037x33=0.825000y33=1.638397y133=1.627882e33=-0.010515x34=0.850000y34=1.654180y134=1.643168e34=-0.011013x35=0.875000y35=1.669842y135=1.658312e35=-0.011530x36=0.900000y36=1.685388y136=1.673320e36=-0.012068x37=0.925000y37=1.700823y137=1.688194e37=-0.012629x38=0.9

36、50000y38=1.716151y138=1.702939e38=-0.013212x39=0.975000y39=1.731376y139=1.717556e39=-0.013820x40=1.000000y40=1.746504y140=1.732051e40=-0.014453thisisimproveEulermethodpleaseinputabandy0:011pleaseinputN:20x0=0.000000y0=1.000000y10=1.000000e0=0.000000x1=0.050000y1=1.048869y11=1.048809e1=-0.000060x2=0.

37、100000y2=1.095561y12=1.095445e2=-0.000116x3=0.150000y3=1.140345y13=1.140175e3=-0.000169x4=0.200000y4=1.183437y14=1.183216e4=-0.000221x5=0.250000y5=1.225017y15=1.224745e5=-0.000272x6=0.300000y6=1.265236y16=1.264911e6=-0.000325x7=0.350000y7=1.304219y17=1.303841e7=-0.000378x8=0.400000y8=1.342074y18=1.3

38、41641e8=-0.000433x9=0.450000y9=1.378896y19=1.378405e9=-0.000492x10=0.500000y10=1.414766y110=1.414214e10=-0.000553x11=0.550000y11=1.449756y111=1.449138e11=-0.000618x12=0.600000y12=1.483927y112=1.483240e12=-0.000687x13=0.650000y13=1.517337y113=1.516575e13=-0.000762x14=0.700000y14=1.550035y114=1.549193

39、e14=-0.000841x15=0.750000y15=1.582066y115=1.581139e15=-0.000927x16=0.800000y16=1.613472y116=1.612452e16=-0.001020x17=0.850000y17=1.644289y117=1.643168e17=-0.001121x18=0.900000y18=1.674551y118=1.673320e18=-0.001230x19=0.950000y19=1.704288y119=1.702939e19=-0.001349x20=1.000000y20=1.733529y120=1.732051

40、e20=-0.0014784.1.4 结论这道题用普通的微分方程也能列式求解，关键在于如何预测若干年后的净资产趋势，用普通的微分方程就无法进行预测，且人工计算量相当大，这里我使用欧拉方法可以计算出精度以及误差，通过电脑运行程序就可以预测出若干年后的净资产情况，欧拉方法的使用使得解题更加方便且精确。4.2 龙格一库塔解题及其数学模型4.2.1 问题提出例4.2.1两种果树寄生虫，其数量分别是u=u(t),v=v(t),其中一种寄生虫以吃另一种寄生虫为生，两种寄生虫的增长函数如下列常微分方程组所示：du于=0.09u(1-)-0.45uv20dv=0.06v(1一看)-0.001uvu(0)=1.

41、6v(0)=1.2问题：预测3年后这一对寄生虫的数量4.2.2 模型建立分析：这是一个典型的常微分方程例题，要求预测出3年后这对寄生虫的数量。确定变量：假定时间是寄生虫数量的积分函数，不妨设变量t为时间(以年为单位)。du=0.09u(1-)-0.45uv由题知有两种寄生虫u和v,u寄生虫年增长函数为dt20,丫寄dv-0.06v(1-)-0.001uv生虫年增长函数为dt15,初始值：u寄生虫为1.6,v寄生虫为1.2,由于其中一种寄生虫以吃另一种寄生虫为生，我们可建立u和v的关联函数f(u,v),g(u,v)。归纳后得到公式：f(u,v)=0.09u(1-)-0.45uv20(4-2),、

42、v、g(u,v)=0.06v(1-一)-0.001uvL15(4-2)即为例4.2.1所述问题的微分方程模型。Inf(un,vn)+hy/【g(un,vn)j4.2.3 解决问题在本题中f(t,u,v)=f(u,v),g(t,u,v)=g(u,v),用Euler预估一校正公式Qnhf(un,vn)-iyn)21f(un,vn)f(un1,vn1)(g(un+l,vn+1)J取h=1,计算结果如表2.表2t/年u(t)v(t)11.61.221.024571.2683430.6409121.336640.3912111.41077我把龙格一库塔的C程序进行编辑后得到结果：(ei表示误差)plea

43、seinputabandy0:：011pleaseinputN:10x0=0.000000y0=1.000000y10=1.000000e0=0.000000x1=0.100000y1=1.095446y11=1.095445e1=-0.000000x2=0.200000y2=1.183217y12=1.183216e2=-0.000001x3=0.300000y3=1.264912y13=1.264911e3=-0.000001x4=0.400000y4=1.341642y14=1.341641e4=-0.000001x5=0.500000y5=1.414215y15=1.414214e5=

44、-0.000002x6=0.600000y6=1.483242y16=1.483240e6=-0.000002x7=0.700000y7=1.549196y17=1.549193e7=-0.000003x8=0.800000y8=1.612455y18=1.612452e8=-0.000004x9=0.900000y9=1.673324y19=1.673320e9=-0.000004x10=1.000000y10=1.732056y110=1.732051e10=-0.0000054.2.4 结论从上面的结果可以分析，用每一种方法计算节点x=0.1、0.20.9、1.0,上的值y都需要计算4次

45、f(x,y)的值，即它们的计算量基本相同，其结果是经典的龙格一库塔方法的精度最好，龙格一库塔方法的推导是基于Taylor级数方法的，因而在使用高阶龙格一库塔方法计算时可以很精确的推算出寄生虫每一年的数量。4.3 分别使用二阶、三阶龙格一库塔法解初值问题对一些特殊的微分方程，使用欧拉方法和低阶的龙格一库塔方法也能达到很高的精度，例如：微分方程的解析解是一次函数，则用欧拉方法求得的数值解与准确解相符，微分方程的解析解是二次函数，则用二阶龙格一库塔方法求得的数值解与准确解相符。微分方程的解析解是三次多项式，用三阶龙格一库塔方法求得的数值解与准确解相符。4.3.1 建立模型建立初值问题模型:Mx1y(

46、0)=1(0x1)(4.1)4.3.2 用不同的龙格一库塔法解(4.1)初值问题模型(1)用二阶龙格库塔方法求解初值问题(4.1)结果如下：thisisthesecond-orderRunge-Kutta(Heun)methodpleaseinputabandy0:011pleaseinputN:10x0=0.000000y0=1.000000y10=1.000000e0=0.000000x1=0.100000y1=1.110000y11=1.110000e1=0.000000x2=0.200000y2=1.240000y12=1.240000e2=0.000000x3=0.300000y3=

47、1.390000y13=1.390000e3=0.000000x4=0.400000y4=1.560000y14=1.560000e4=0.000000x5=0.500000y5=1.750000y15=1.750000e5=0.000000x6=0.600000y6=1.960000y16=1.960000e6=0.000000x7=0.700000y7=2.190000y17=2.190000e7=0.000000x8=0.800000y8=2.440000y18=2.440000e8=0.000000x9=0.900000y9=2.710000y19=2.710000e9=0.00000

48、0x10=1.000000y10=3.000000y110=3.000000e10=0.000000thisisimproveEulermethodpleaseinputabandy0：011pleaseinputN:10x0=0.000000y0=1.000000y10=1.000000e0=0.000000x1=0.100000y1=1.110000y11=1.110000e1=0.000000x2=0.200000y2=1.240000y12=1.240000e2=0.000000x3=0.300000y3=1.390000y13=1.390000e3=0.000000x4=0.4000

49、00y4=1.560000y14=1.560000e4=0.000000x5=0.500000y5=1.750000y15=1.750000e5=0.000000x6=0.600000y6=1.960000y16=1.960000e6=0.000000x7=0.700000y7=2.190000y17=2.190000e7=0.000000x8=0.800000y8=2.440000y18=2.440000e8=0.000000x9=0.900000y9=2.710000y19=2.710000e9=0.000000x10=1.000000y10=3.000000y110=3.000000e1

50、0=0.000000thisisthesecond-orderRunge-Kutta(middle)methodpleaseinputabandy0：011pleaseinputN:10x0=0.000000y0=1.000000y10=1.000000e0=0.000000x1=0.100000y1=1.110000y11=1.110000e1=0.000000x2=0.200000y2=1.240000y12=1.240000e2=0.000000x3=0.300000y3=1.390000y13=1.390000e3=0.000000x4=0.400000y4=1.560000y14=1

51、.560000e4=0.000000x5=0.500000y5=1.750000y15=1.750000e5=0.000000x6=0.600000y6=1.960000y16=1.960000e6=0.000000x7=0.700000y7=2.190000y17=2.190000e7=0.000000x8=0.800000y8=2.440000y18=2.440000e8=0.000000x9=0.900000y9=2.710000y19=2.710000e9=0.000000x10=1.000000y10=3.000000y110=3.000000e10=0.000000(2)用三阶龙格

52、一库塔方法求解初值问题模型:-dy=3x2+2x+1dxy(0)=1(0x三2)(4.2)结果如下：thisisthethird-orderHeunmethod:pleaseinputabandy0:021pleaseinputN:10x0=0.000000y0=1.000000y10=1.000000e0=0.000000x1=0.200000y1=1.248000y11=1.248000e1=0.000000x2=0.400000y2=1.624000y12=1.624000e2=0.000000x3=0.600000y3=2.176000y13=2.176000e3=0.000000x4=0.800000y4=2.952000y14=2.952000e4=0.000000x5=1.000000y5=4.000000y15=4.000000e5=-0.000000x6=1.200000y6=5.368000y16=5.36800