数模微分方程的数值求解技术

1、常微分方程的数值求解技术1 Euler方法2 Runge-Kutta 方法3线性多步方法1 Euler方法1.1 从一个人口增长(Malthus )模型谈起英国人口统计学家马尔萨斯(Malthus,1766-1834年)调查了英国100多年人口出生统计资料，发现人口出生率稳定于一个常数。并提出了著名的Malthus人口增长模型。此模型的基本假设是：在考虑一个国家或地区的人口总数随时间变化的人口自然增长过程中，略去迁移、 自然环境条件等因素对人口变化的影响，视净相对增长率(出生率与死亡率之差)是常数。 即单位时间内人口的增长量与人口成正比，比例系数为r。此假设就是一个常微分方程的问题。记x(t)

2、为时刻t的人口，视x(t)为连续、可微函数。据 Malthus的假设，在t到t ：二t时间内人口的增长量为 x(t Ct) -x(t)二rx(t)，又设t =t0时的人口为x0，令4 - 0时,可得x(t)满足方程匸dxrxdtx(0) = x假设取1790年为初始状态t0 =0，已知X。=3.9 10与 1800年的实际人口6x(1) = 5.3 10，由最小二乘法拟合， 可以定出r =0.307。计算1790年1880年间的人口 (百 万)，h = 1。即求解微分方程=0.307xx(0) =3.9以下先介绍Euler方法，再给出此微分方程的数值解。1). Euler 公式yn 1 =yn

3、hf (xn, yn) (n 二 0,1,2, )(2)2).改进的Euler公式yn 1 二 Yn hf(xn,Yn)y=dsolve( Dy=0.307*y , y(0)=3.9 , x)运行后显示精确解y=39/10*exp(307/1000*x);(1) 建立并保存名为fun 11.m的M文件函数。fun cti on f=fun 1(x,y)f=0.307*y;(2) 建立并保存名为 Euler_f.m 和Euler_r.m的M文件函数。(3) 在Matlab工作窗口输入程序% 求精确解 y= dsolve( Dy=0.307*y , y(0)=3.9 , x)% y=39/10*e

4、xp(307/1000*x);% y1为用欧拉公式计算的数值解；% y2为用改进的欧拉公式计算的数值解.x=0:1:9;y=39/10*exp(307/1000*x);x,y1=Euler_f( fun 1,0,3.9,1,9)x,y2=Euler_r( fun 1,0,3.9,1,9)plot(x,y1,*r,x,y2,pb,x,y)legend(用欧拉公式计算 dy/dx=0.307y,y(0)=3.9 在0,9上的数值解,用改进的欧拉公式计算 dy/dx=0.307y,y(0)=3.9 在0,9上的数值解,dy/dx=0.307y,y(0)=3.9 在0,9 上的精确解)运行后的计算结果

5、如下表6-2表6-2用Euler公式及改进的Euler公式求17901880年的人口数据(10 6)年份1790180018101820183018401850186018701880Xi0123456789屮3.90005.09736.66228.707511.380614.874519.441025.409433.210043.4055y23.90005.28117.15129.683713.112917.756524.044532.559344.089359.7024y3.90005.30147.20659.796013.316118.101224.605733.447545.46656

6、1.8045在计算机中显示其数值解与精确解的图形如下图6-2图6-2用Euler公式和改进的Euler公式求dy/dx=0.307y,y(0)=3.9 在0,9上的数值解与精确解图形例3 用Euler公式和改进的Euler公式求初值问题_3y+8x_7(0 兰 x 兰1)J(0) =1的数值解(取h =0.1)，并作出此数值解和该微分方程的精确解的图像。 解具体的计算公式为yi 比=yj + 0.1(3yj + 8xi - 7) = 0.7y 0.8xi - 0.7 Euler公式lyo =1改进的Euler公式= 0.7%0.8备-0.70.1 _二 y2【3yi 8Xi -0.7 (-3N

7、 1 8x 1 -0.7)=y 0.1 -3yi 8为-0.7 (-3(0.7%0.8xi-0.7) 8xi 0.7)= 0.745% 0.28x0.4xi d 0.035y 二1利用Matlab编程过程如下: 求精确解y=dsolve(Dy=8*x-3*y-7,y(0)=1,x)运行后显示精确解y=8/3*x-29/9+38/9*exp(-3*x)(1) 建立并保存名为fun 12.m的M文件函数。function f=fun 12(x,y)f=8*x-3*y-7;(2) 建立并保存名为 Euler_f.m 和Euler_r.m的M文件函数。(3) 在Matlab工作窗口输入程序% 求精确解

8、 y= dsolve( Dy=-3*y+8*x-7 , (0)=1 , x)% y=8/3*x-29/9+38/9*exp(-3*x)为精确解;% y1为Euler公式计算的数值解；% y2为改进的Euler公式计算的数值解.x=0:0.1:1;y=8/3*x-29/9+38/9*exp(-3*x)x,y1=Euler_f( fun 12,0,1,0.1,10)x,y2=Euler_r( fun 12,0,1,0.1,10)plot(x,y1,*r,x,y2,pb,x,y)legend(用欧拉公式计算 dy/dx=-3y+8x-7,y(0)=1 在0,1上的数值解,用改进的欧拉 公式计算 dy

9、/dx=-3y+8x-7,y(0)=1 在0,1上的数值解,dy/dx=-3y+8x-7,y(0)=1 在 0,1上的精确解)运行后的计算结果如下表6-3表6-3 用Euler公式及改进的Euler公式求初值问题dy/dx=-3y+8x-7,y(0)=1 在0,1上的数值解及精确解Xny1y2y01.00001.00001.00000.100000.19000.17230.2000-0.6200-0.3455-0.37170.3000-0.9740-0.6764-0.70560.4000-1.1418-0.8549-0.88380.5000-1.1793-0.9199-0.94680.6000

10、-1.1255-0.9003-0.92430.7000-1.0078-0.8177-0.83850.8000-0.8455-0.6882-0.70590.9000-0.6518-0.5237-0.53851.0000-0.4363-0.3332-0.3453在计算机中显示其精确解和数值解的图形如下图6-3+月欧施鱼式计算dy/dx=-3y舷7，沌尸1祜0.1 上苗数咕堀ft用比进的陆摄仑或计算岬d疔谢钊=1在|0,11的数值削 dy/d沪-刊匚在上的:瘠确解D.5A1D-0 5 I1 + * +iiii1i1hi0.10.20.30.40.50.G0.70.B0.9图6-3用Euler公式及改

11、进的Euler公式求dy/dx=-3y+8x-7,y(0)=1 在0,1上的数值解与精确解图形 由例2及例3的数据及数值解与精确解的图形可以看出，改进的Euler公式比Euler公式的精度要高。2 Run ge-Kutta 方法2.1 从一个阻滞增长(Logistic )模型谈起从上面提出的人口指数增长(Malthus )模型可知，作短期人口预测可以得到较好的结果，但由模型的精确解 x(t) =X0e可知，lim x(t) = lim x0ert = :，即用指数模型作长期预测人 t0口，人口将无限增长，这是不合理的。事实上，自然资源、环境条件等因素对人口的增长起 着阻滞作用，并且随着人口的增

12、加，阻滞作用会越来越大，于是对指数增长模型进行一种改 进得到阻滞增长(Logistic )模型。荷兰生物数学家Verhulst引入常数xm，用来表示自然资源和环境条件所能容许的最大人 口数，并假定净增长率为r(x(t)=r(1 一晋)即净增长率随着x(t)的增加而减小。当x(t)r xm时，净增长率趋于零。其中r, xm是根据人口统计数据或经验确定的常数，比例因子1_也体现了对人口增长的阻滞作用。上式还可解Xm释为净增长率r(x(t)与人口尚未实现部分 (对最大容许量而言)的比例xm ft)成正比，比例系数为固有增长率r(r 0)。所以由假设，指数增长模型被改为上式通常被称为阻滞增长（dx(t

13、)dt二 rx (t)(1也)XmLogistic ）模型。写为带初值问题的微分方程为Idx 一 x、rx（1）“ dtXm公（0）= X。假设取1790年为初始状态6 6t0 =0，已知 x。=3.9 10，取 r = 0.3134，xm =197 10 。计算1790年1880年间的人口（ 106），即求解微分方程齐a3%话)xx(0) =3.9其中0 ex乞9, h=1。1. Runge-Kutta 格式.四级四阶R-K方法在四级R-K方法中，最常用的是标准的（或经典的）四级四阶方法yn 1 Tn (& 2k2 2k3 k4)k1 = f (XnWn)11二 f (xnh, ynhk12

14、211二 f(xnh, ynhk2)22f(xn h, yn hk3k2k3k4(6.3.10亦可写成如下等价形式yn-1k1k2k3k41=yn -(k1 2k2 2k3 k4) hf (xn, yn)112h, yn -k1)11尹2)yn k3=hf-hf=hf(xn(xn(xn尹ynh,(6.3.112 .阻滞增长（Logistic 例1 以上提出的）模型的题解阻滞增长(Logistic )模求解微分方程dx .3134（）y，试用经典的四阶R-K公式计算此微分方程的数值解。其中，1790年,y（0） =3.91880年转化为自变量取值 0岂x乞9，步长h =1。解 由公式(6.3.1

15、0 )，此问题的经典的四阶 R-K公式的具体格式为x=x +ih, i =0,1,2,川,9,0.3134 2ki = f (Xn, yn) =0.3134yn -一yn,197k2 = f (xn +0.5h,yn +0佛&) =0.3134( yn “册)03134“。*2, 197“0.31342k3 = f(xn +0.5h,yn +0.5hk2)=0.3134(yn +0.5k2)(yn+0.5k2)2,1970.31342kq = f 区 +h, yn +hk3)=0.3134(yn +kj -Wn*),197hyn+ = y-(k1 +2kH2k y=dsolve( Dy=0.3

16、134*y*(1-y/197), y(0)=3.9,x)运行后输出结果为精确解y如下y=197心+1931/39*exp(-1567/5000*x)利用Matlab编程过程如下：(1) 建立并保存名为fun31.m的M文件函数。function f=fun 31(x,y)f=0.3134*(1-1/197*y).*y;(2) 建立并保存名为 Runge_Kutta4.m的M文件函数。(3) 在Matlab工作窗口输入程序% y=dsolve( Dy=0.3134*y*(1-1/197*y), y(0)=3.9,x);% y=197/(1 + 1931/39*exp(-1567/5000*x)

17、为精确解；% y1为用经典的四阶R-K公式计算的数值解.x=0:1:9;y=197./(1+1931./39*exp(-1567./5000*x)x,y1=R un ge_Kutta4( fun 31,0,3.9,1,9)plot(x,y1,pb, x,y ，-r)legend(用经典的四阶 R-K公式求 dy/dx=0.3134(1-y/197),y(0)=3.9在0,9上的数值解,dy/dx=0.3134(1-y/197),y(0)=3.9在0,9上的精确解)运行后计算结果如下表6-6表6-6用四阶R-K方法估计17901880年份的人口数据(106)年份17901800181018201

18、83018401850186018701880x0123456789y13.90005.29687.17529.686213.015717.383423.033330.210839.121949.8756y3.90005.29697.17559.686713.016517.384623.035230.213339.125349.8799其数值解和精确解的图形如下图6-4图6-4 用四阶R-K方法计算dy/dx=0.3134(1-y/197),y(0)=3.9 在0,9上的数值解与精确解图形例3 用经典的四阶Runge-Kutta公式求初值问题dx 1x2y(0)=0的数值解(取h =0.25)

19、，并作出此数值解和该微分方程的精确解的图形。解先求精确解。输入程序Y=dsolve( (Dy)+(2*x*y)/(1+xA2)-1=0, y(0)=0 , x)运行后输出结果为精确解y如下y= (x+1/3*xA3”(1+xA2)利用Matlab编程过程如下：(1) 建立并保存名为fun32.m的M文件函数。Fun ction f=fun 32(x,y) f=1-(2*x*y)/(1+xA2);(2) 建立并保存名为 Runge_Kutta4.m的M文件函数。(3) 在Matlab工作窗口输入程序% 求精确解 y=dsolve( (Dy)+(2*x*y)/(1+xA2)-1=0, y(0)=0

20、 , x);% y= (x+1/3*xA3)/(1+xA2);% y1为用经典的四阶R-K公式计算的数值解；x=0:0.25:2;y=(x+1/3*x.A3)./(1+x.A2)x,y1=Runge_Kutta4( fun32,0,0,0.25,8)plot(x,y1,pb, x,y,-r)legend(用经典四阶 R-K公式计算 dy/dx=1-(2*x*y)/(1+xA2),y(0)=0 在0,2上的数值解,dy/dx=1-(2*x*y)/(1+xA2),y(0)=0在0,2的精确解)运行后计算结果如表6-8表6-8用四阶 R-K方法计算 dy/dx=1-(2*x*y”(1+xA2),y(

21、0)=0的结果x00.25000.50000.75001.00001.25001.50001.75002.0000y5.29687.17529.686213.015717.383423.033330.210839.121949.8756yi5.29697.17559.686713.016517.384623.035230.213339.125349.8799其数值解和精确解的图形如下图6-51 111111* .用鮭負四阶艮式老听川“ 1 A神1 +爲丫（0冋在|0 21的数值娇 y=dsolve( DY)+0.9*Y-0.9*1.5=0, Y(0)=1 , x)运行后输出结果为精确解如下y=

22、3/2-1/2*exp(-9/10*x)(1) 建立并保存名为fun41.m的M文件函数。function f=fun 41(x,y)f=0.9*(1.5-y);(2) 建立并保存名为 Adamsyx.m的M文件函数。(3) 在Matlab工作窗口输入程序% y=dsolve(Dy=0.9*(1.5-y),y(0)=1,x);% y =3/2-1/2*exp(-9/10*x);% y1为四阶亚当斯预测-校正公式计算的数值解.x=0:1:11;y=3/2-1/2*exp(-9/10*x)x,y1=Adamsyx( fun41,0,1,1,11)plot(x,y1,pb, x,y,-r)legen

23、d(四阶亚当斯预测-校正公式求dy/dx=0.9*(1.5-y),y(0)=1 在0,11上数值解 精确解 y=3/2-1/2*exp(-9/10*x)运行后计算结果如下表6-12表6-12用四阶Adams!测-校正公式估计112月份的产品价格月份xn123456x012345y1.00001.29671.41741.46641.48631.4944y11.00001.29461.41561.46531.49331.5021月份Xn789101112x67891011y1.49771.49911.50161.50091.49961.5000y11.49981.49971.50161.50091

24、.49961.5000其数值解与精确解的图形如下图6-617*四斷戎会式求旳闭円.宁门.$巧片1 ft|D,11裁值解 帝确解y=3-i g叩円/i CTn)1 .b1.51,41.31.21.1-工一丄一 n s r a总 *-ST / /-/ - !.1 ./iL1111124EB101：图6-6 用四阶 Adams!测-校正公式求 dy/dx=0.9*(1.5-y),y(0)=1 在0,11上的数值解与精确解图形例3用四阶Adams预测-校正公式求初值问题dx 1 x2y(0)=0 x+x3的数值解(取h =0.2)，并作出此数值解和该微分方程的精确解y =3的图形。1 +x2解 (1 )建立并保存名为fun42.m的M文件函数。function f=funfcn4 2(x,y)f=1-(2*x*y)/(1+xA2);(2) 建立并保存名为 Adamsyx.m的 M文件函数。(3) 在Matlab工作窗口输入程序% 精确解 y=dsolve(Dy=