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

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

2、视为连续、可微函数。据Malthus的假设，在到时间内人口的增长量为，又设时的人口为，令时，可得满足方程假设取1790年为初始状态，已知与1800年的实际人口，由最小二乘法拟合，可以定出。计算1790年1880年间的人口（百万），。即求解微分方程以下先介绍Euler方法，再给出此微分方程的数值解。1） Euler公式 （1.2） 2）.改进的Euler公式 (1.5)此公式称为改进的Euler公式。它也可写成一个显式公式 (1.6)或 (1.7)3）人口增长（Malthus）模型的题解例1 取1790年为初始状态，已知与1800年的实际人口，由最小二乘法拟合，可以定出。用改进的Euler公式计

3、算1790年1880年间的人口（百万），即求解微分方程其中，。解 改进的Euler公式具体的计算公式为，计算结果如下表6-1。表6-1 用改进的Euler公式估计17901880年份的人口数据（106）年份1790180018101820183018401850186018701880ti0123456789xi3.95.28117.15129.683713.112917.756524.044532.559344.089359.70241.2 软件介绍在本节里，将通过数值例子和Matlab程序介绍本节提出的Euler公式及改进的Euler公式的上机实现方法。1.根据式（6.1.2）及算法6-1

4、所表示的Euler公式，编写Matlab程序(函数名：Euler_f.m)。function x,y=Euler_f(fun1,x0,y0,h,N)% fun1 为一阶微分方程的函数；% x0,y0为初始条件； % h为区间步长；% N为区间的个数；% x为Xn构成的向量；% y为Yn构成的向量.x=zeros(1,N+1);y=zeros(1,N+1);x(1)=x0;y(1)=y0;for n=1:N x(n+1)=x(n)+h; y(n+1)=y(n)+h*feval(fun1,x(n),y(n);end2.根据式（6.1.7）及算法6-2所表示的改进的Euler公式，编写Matlab程

5、序(函数名：Euler_r.m)。function x,y=Euler_r(fun1,x0,y0,h,N)% fun1 为一阶微分方程的函数；% x0,y0为初始条件； % h为区间步长；% N为区间的个数；% x为Xn构成的向量；% y为Yn构成的向量.x=zeros(1,N+1);y=zeros(1,N+1);x(1)=x0;y(1)=y0;for n=1:N x(n+1)=x(n)+h; ybar=y(n)+h*feval(fun1,x(n),y(n); y(n+1)=y(n)+h/2*(feval(fun1,x(n),y(n)+feval(fun1,x(n+1),ybar);end例2

6、 采用Euler公式（Euler_f.m）和改进的Euler公式（Euler_r.m）的Matlab程序上机实现人口指数增长（Malthus）模型所反映的微分方程的题解。其中，。解 利用Matlsb编程过程如下：求精确解>>y=dsolve(Dy=0.307*y,y(0)=3.9,x)运行后显示精确解y=39/10*exp(307/1000*x);（1）建立并保存名为fun11.m的M文件函数。function f=fun1(x,y)f=0.307*y;（2）建立并保存名为Euler_f.m和Euler_r.m的M文件函数。（3）在Matlab工作窗口输入程序% 求精确解y= ds

7、olve(Dy=0.307*y,y(0)=3.9,x)% y=39/10*exp(307/1000*x);% y1为用欧拉公式计算的数值解;% y2为用改进的欧拉公式计算的数值解.x=0:1:9;y=39/10*exp(307/1000*x); x,y1=Euler_f( fun1,0,3.9,1,9) x,y2=Euler_r( fun1,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上的数值解','用改进的欧拉公式计算d

8、y/dx=0.307y,y(0)=3.9在0,9上的数值解','dy/dx=0.307y,y(0)=3.9在0,9上的精确解')运行后的计算结果如下表6-2表6-2 用Euler公式及改进的Euler公式求17901880年的人口数据 (106)年份1790180018101820183018401850186018701880xi0123456789y13.90005.09736.66228.707511.380614.874519.441025.409433.210043.4055y23.90005.28117.15129.683713.112917.756524.

9、044532.559344.0893 59.7024y3.90005.30147.20659.796013.316118.101224.605733.447545.466561.8045在计算机中显示其数值解与精确解的图形如下图6-2图6-2 用Euler公式和改进的Euler公式求dy/dx=0.307y,y(0)=3.9在0,9上的数值解与精确解图形例3 用Euler公式和改进的Euler公式求初值问题的数值解（取），并作出此数值解和该微分方程的精确解的图像。解 具体的计算公式为Euler公式改进的Euler公式利用Matlab编程过程如下：求精确解>>y=dsolve(

10、9;Dy=8*x-3*y-7','y(0)=1','x')运行后显示精确解y=8/3*x-29/9+38/9*exp(-3*x)（1）建立并保存名为fun12.m的M文件函数。function f=fun12(x,y)f=8*x-3*y-7;（2）建立并保存名为Euler_f.m和Euler_r.m的M文件函数。（3）在Matlab工作窗口输入程序% 求精确解y= dsolve(Dy=-3*y+8*x-7,y(0)=1,x)% y=8/3*x-29/9+38/9*exp(-3*x)为精确解;% y1为Euler公式计算的数值解;% y2为改进的Euler

11、公式计算的数值解.x=0:0.1:1;y=8/3*x-29/9+38/9*exp(-3*x) x,y1=Euler_f( fun12,0,1,0.1,10) x,y2=Euler_r( fun12,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/dx=-3y+8x-7,y(0)=1在0,1上的数值解','dy/dx=-3y+8x-7,y(0)=1在0,1上的精确解&#

12、39;)运行后的计算结果如下表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-1.1255-0.9003-0.92430.7000-1.0078-0.8177-0.83850.8000-0.84

13、55-0.6882-0.70590.9000-0.6518-0.5237-0.53851.0000-0.4363-0.3332-0.3453在计算机中显示其精确解和数值解的图形如下图6-3图6-3 用Euler公式及改进的Euler公式求dy/dx=-3y+8x-7,y(0)=1在0,1上的数值解与精确解图形由例2及例3的数据及数值解与精确解的图形可以看出，改进的Euler公式比Euler公式的精度要高。2 Runge-Kutta方法2.1 从一个阻滞增长（Logistic）模型谈起从上面提出的人口指数增长（Malthus）模型可知，作短期人口预测可以得到较好的结果，但由模型的精确解可知，即用

14、指数模型作长期预测人口，人口将无限增长，这是不合理的。事实上，自然资源、环境条件等因素对人口的增长起着阻滞作用，并且随着人口的增加，阻滞作用会越来越大，于是对指数增长模型进行一种改进得到阻滞增长（Logistic）模型。荷兰生物数学家Verhulst引入常数，用来表示自然资源和环境条件所能容许的最大人口数，并假定净增长率为即净增长率随着的增加而减小。当时，净增长率趋于零。其中是根据人口统计数据或经验确定的常数，比例因子体现了对人口增长的阻滞作用。上式还可解释为净增长率与人口尚未实现部分（对最大容许量而言）的比例成正比，比例系数为固有增长率。所以由假设，指数增长模型被改为上式通常被称为阻滞增长（

15、Logistic）模型。写为带初值问题的微分方程为假设取1790年为初始状态，已知，取，。计算1790年1880年间的人口（106），即求解微分方程其中，。1 Runge-Kutta格式. 四级四阶R-K方法在四级R-K方法中，最常用的是标准的（或经典的）四级四阶方法 （6.3.10）亦可写成如下等价形式 （6.3.11）2阻滞增长（Logistic）模型的题解例1 以上提出的阻滞增长（Logistic）模型转化为求解微分方程，试用经典的四阶R-K公式计算此微分方程的数值解。其中， 1790年1880年转化为自变量取值，步长。解 由公式（6.3.10），此问题的经典的四阶R-K公式的具体格式为

16、具体的数值解见下表6-5表6-5 用四阶R-K方法估计17901880年份的人口数据（106）年份1790180018101820183018401850186018701880xn0123456789yn3.90005.29687.17529.686213.015717.383423.033330.210839.121949.87562.2 软件介绍在本节里，将通过数值例子和Matlab程序介绍本节提出的经典四阶Runge-Kutta公式的上机实现方法。根据式（6.3.10）及算法6-3所表示的经典四阶Runge-Kutta公式，编写Matlab程序(函数名：Runge_Kutta4.m)。

17、function x,y=Runge_Kutta4(fun,x0,y0,h,N)% fun 为一阶微分方程的函数；% x0,y0为初始条件； % h为区间步长；% N为区间的个数；% x为Xn构成的向量；% y为Yn构成的向量；x=zeros(1,N+1);y=zeros(1,N+1);x(1)=x0;y(1)=y0;for n=1:N x(n+1)=x(n)+h; k1=h*feval(fun,x(n),y(n); k2=h*feval(fun,x(n)+1/2*h,y(n)+1/2*k1); k3=h*feval(fun,x(n)+1/2*h,y(n)+1/2*k2); k4=h*feva

18、l(fun,x(n)+h,y(n)+k3); y(n+1)=y(n)+1/6*(k1+2*k2+2*k3+k4);end例2 用经典的四阶Runge-Kutta公式的Matlab程序上机实现阻滞增长（Logistic）模型的数值解，并作出此数值解和该微分方程的精确解的图形。解 先求精确解。输入程序>> y=dsolve(Dy=0.3134*y*(1-y/197)',y(0)=3.9',x')运行后输出结果为精确解y如下y= 197/(1+1931/39*exp(-1567/5000*x)利用Matlab编程过程如下：（1）建立并保存名为fun31.m的M文件

19、函数。function f=fun31(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)为精确解;% y1为用经典的四阶R-K公式计算的数值解.x=0:1:9;y=197./(1+1931./39*exp(-1567./5000*x)x,y1=Runge_Kutta4( fun31,0,

20、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）年份1790180018101820183018401850186018701880x0123456789y13.90005.29687.17529.686213.015717.

21、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公式求初值问题的数值解（取），并作出此数值解和该微分方程的精确解的图形。解 先求精确解。输入程序>>Y=dsolve(Dy)+(2*x*y)/(1+x2)-1=0,y(0)=0,x)运行后输出结果

22、为精确解y如下y= (x+1/3*x3)/(1+x2)利用Matlab编程过程如下：（1）建立并保存名为fun32.m的M文件函数。Function f=fun32(x,y)f=1-(2*x*y)/(1+x2);（2）建立并保存名为Runge_Kutta4.m的M文件函数。（3）在Matlab工作窗口输入程序% 求精确解y=dsolve(Dy)+(2*x*y)/(1+x2)-1=0,y(0)=0,x);% y= (x+1/3*x3)/(1+x2);% y1为用经典的四阶R-K公式计算的数值解;x=0:0.25:2;y=(x+1/3*x.3)./(1+x.2)x,y1=Runge_Kutta4(

23、 fun32,0,0,0.25,8)plot(x,y1,'pb', x,y,'-r')legend('用经典四阶R-K公式计算dy/dx=1-(2*x*y)/(1+x2),y(0)=0在0,2上的数值解','dy/dx=1-(2*x*y)/(1+x2),y(0)=0在0,2的精确解')运行后计算结果如表6-8表6-8 用四阶R-K方法计算dy/dx=1-(2*x*y)/(1+x2),y(0)=0的结果x00.25000.50000.75001.00001.25001.50001.75002.0000y5.29687.17529.6

24、86213.015717.383423.033330.210839.121949.8756y15.29697.17559.686713.016517.384623.035230.213339.125349.8799其数值解和精确解的图形如下图6-5图6-5 用经典的四阶R-K公式计算dy/dx=1-(2*x*y)/(1+x2),y(0)=0在0,2上的数值解与数值解图形由例2及例3的数据及数值解与精确解的图形可以看出，经典的四阶R-K公式的精度较高，数值解几乎与精确解重合。3 线性多步法3.1 从一个市场动态均衡价格模型谈起设某种商品的需求量与供给量只与该商品的价格有关。可设其需求函数与供给函

25、数分别为则该商品的均衡价格为。一般地，当市场上该商品供过于求（）时，价格将下跌；当供不应求（）时，价格将上涨。因此，该商品在市场上的价格将随着时间的变化而围绕着均衡价格上下波动，价格是时间的函数。根据上述供求关系变化影响价格变化的分析，可以假设时刻价格的变化率与时刻的超额需求量成正比，即设其中为正的常数，用来反映价格的调整速度。将需求函数与供给函数代入上式，整理得市场动态均衡价格模型其中。假设需求函数，供给函数，反映价格的调整速度的常数。则均衡价格，常数。所以市场动态均衡价格模型为以下先介绍线性多步法的Adams（阿当姆斯）方法，再来求解市场动态均衡价格模型所反映的此微分方程的数值解。1）Ad

26、ams（阿当姆斯）方法的显式与隐式公式形如 （6.4.7）的步法，称为Adams方法。时为显式方法；时为隐式方法。 表6-9 Adams显式公式kp公式12341234表6-10 Adams隐式公式kp公式12342345注：为表现出方法的阶数，以下记表6-10中的为。最常用的是四阶Adams显式公式和四阶Adams隐式公式。2） Adams预测-校正方法 此公式称为四阶的Adams预测-校正公式。为了编程时书写方便，四阶的Adams预测-校正公式常记为 （6.4.10）因为式（6.4.10）是四阶公式，所以它的起步值除了给定的以外，通常用经典的四阶R-K公式计算。3） 市场动态均衡价格模型的

27、题解例1 若市场动态均衡价格模型记为假设一年中第一个月的为初始状态且，试用四阶Adams预测-校正公式求一年中每个月的价格情况。其中，。解 由四阶的Adams预测-校正公式（6.4.10）知，迭代的初始值由经典的四阶R-K方法计算，由四阶的Adams预测-校正公式计算。由公式（6.3.10），此问题的经典的四阶R-K公式的具体格式为由公式（6.4.10），此问题的四阶Adams预测-校正公式的具体格式为具体的数值解见下表6-11表6-11 用四阶Adams预测-校正公式估计112月份的产品价格月份xn123456yn1.00001.29461.41561.46531.49331.5021月份x

28、n789101112yn1.49981.49971.50161.50091.49961.50003.2 软件介绍在本节里，将通过数值例子和Matlab程序介绍本节提出的四阶Adams预测-校正公式的上机实现方法。根据式（6.4.10）及算法6-5所表示的四阶Adams预测-校正公式，编写Matlab程序(函数名：Adamsyx.m)。function x,y=Adamsyx(fun4,x0,y0,h,N)% fun4为一阶微分方程的函数；% x0,y0为初始条件； % h为区间步长；% N为区间的个数；% x为Xn构成的向量；% y为Yn构成的向量.x=zeros(1,N+1);y=zeros

29、(1,N+1);x(1)=x0;y(1)=y0;for n=1:N x(n+1)=x(n)+h; if n<4 k1=h*feval(fun4,x(n),y(n); k2=h*feval(fun4,x(n)+1/2*h,y(n)+1/2*k1); k3=h*feval(fun4,x(n)+1/2*h,y(n)+1/2*k2); k4=h*feval(fun4,x(n)+h,y(n)+k3); y(n+1)=y(n)+1/6*(k1+2*k2+2*k3+k4); else k1=feval(fun4,x(n),y(n); k2=feval(fun4,x(n-1),y(n-1); k3=fe

30、val(fun4,x(n-2),y(n-2); k4=feval(fun4,x(n-3),y(n-3); ybar=y(n)+h/24*(55*k1-59*k2+37*k3-9*k4); k0=feval(fun4,x(n+1),ybar); y(n+1)=y(n)+h/24*(9*k0+19*k1-5*k2+k3);endend例2 试用四阶Adams预测-校正公式编写的Matlab程序上机实现市场动态均衡价格模型的题解。其中。即估计一年中112月份的价格情况。解 先求精确解。输入程序>>y=dsolve(DY)+0.9*Y-0.9*1.5=0,Y(0)=1,x)运行后输出结果为

31、精确解如下y=3/2-1/2*exp(-9/10*x)（1）建立并保存名为fun41.m的M文件函数。function f=fun41(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,

32、0,1,1,11)plot(x,y1,'pb', x,y,'-r')legend('四阶亚当斯预测-校正公式求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.