微分方程数值解

微

分方程数值解及其应用

绪论 自然界中的许多事物的运动和变化规律都可以用微分方程来描述，因此对工程和科学技术中的实际问题的研究中, 常常需要求解微分方程．但往往只有少数较简单和典型的微分方程可求出其解析解，在大多数情况下,只能用近似法求解，数值解法是一类重要的近似方法．本文主要讨论一阶常微分方程的初值问题的数值解法，探讨这些算法在处理来自生活实际问题中的应用，并结合MATLAB 软件，动手编程予以解决． １ 微分方程的初值问题[1]

1.1 预备知识

在对生活实际问题的研究中，通常需要考虑一阶微分方程的初值问题

00(,)()dy f x y dx y x y ?=???=? （1）

这里(),f x y 是矩形区域R ：00,x x a y y b -≤-≤上的连续函数．

对初值问题（1）需要考虑以下问题：方程是否一定有解呢?若有解，有多少个解呢?下面给出相关的概念与定理．

定义1 Lipschitz 条件[1][2]：矩形区域R ：00,x x a y y b -≤-≤上的连续函数(),f x y 若满足：存在常数0L >，使得不等式()()1212,,f x y f x y L y y -≤-对所有()()12,,,x y x y R ∈都成立，则称(),f x y 在R 上关于y 满足Lipschitz 条件. 定理 1 解的存在唯一性定理[1][3]：设f 在区域()}{,,D x y a x b y R =≤≤∈上连续，关于y 满足Lipschitz 条件，则对任意的[]00,,∈∈x a b y R ，常微分方程初值问题（1）当[],x a b ∈时存在唯一的连续解()y x ．

该定理保证若一个函数(),f x y 关于y 满足Lipschitz 条件，它所对应的微分方程的初值问题就有唯一解.在解的存在唯一性得到保证的前提下，自然要考虑方程的求

解问题．求解微分方程虽然有多种解析方法，但根据工程和科学实践问题所得到的微分方程往往很复杂，在很多情况下不能或很难给出解析解，有时即使能求出形式解，也往往因形式过于复杂或计算量太大而不实用，因此从实际问题中归结出来的微分方程主要依靠数值解法．

定义 2 微分方程数值解：对初值问题（1）寻求数值解就是寻求解()y x 在一系列离散节点

上的近似解0121,,,,,,n n y y y y y +L L ，相邻两个节点的间距1n n n h x x +=-称为步长.在一般情况下假定()0,1,i h h i ==L 为常数，这时节点为0,0,1,2,n x x nh n =+=L ．

要求微分方程数值解，首先要建立数值算法，即对初值问题（1）中的方程离散化，建立求解数值解法的递推公式．一类是计算1n y +时只用到前一点的值n y ，称为单步法；另一类是用到1n y +前面k 点的值11,,,n n n k y y y --+L 称为k 步法.

对初值问题（1）式的单步法可用一般形式表示为

11(,,,)n n n n n y y h x y y h ?++=+?，

其中多元函数?与(),f x y 有关，当?含有1n y +时，方法是隐式的；若?中不含1n y +，则为显式方法，所以显式单步法可表示为

1(,,)n n n n y y h x y h ?+=+?.

（2）

设()y x 是初值问题（1）的准确解，称()()()11(,,)n n n n n T y x y x h x y x h ?++=--?为显式单步法（2）的局部截断误差. 若存在最大正整数p ,使显式单步法（2）式的局部截断误差满足()()()()11,,p n T y x h y x h x y h O h ?++=+--=，则称（2）式有p 阶精度.

1.2几种常用的数值解法及其分析、比较

1.2.1欧拉法与后退欧拉法

1）欧拉法：欧拉曾简单地用差分代替微分，即利用公式

将初值问题（1）离散化，则问题（1）可化为

1(,),n n n n y y h f x y +=+?0n x x n h =+?， （3）

此方法称为欧拉法.

欧拉方法的几何意义在数值计算公式中体现了出来.在xy 平面上，一阶微分方程的解()y y x =称作它的积分曲线.积分曲线上一点(),x y 的切线斜率等于函数(),f x y ，按函数(),f x y 在xy 平面上建立一个方向场，那么，积分曲线上每一点的切线方向均与方向场在该点的方向相一致.

基于上述几何解释，从初始点000(,)P x y 出发，先依方向场在该点的方向上推进到1x x =上一点1P ，再从1P 依方向场的方向推进到2x x =上一点2P ，循环前进便作出一条折线012P PP L ，因此欧拉方法又称为折线法.若初值0y 已知，则由（3）式可逐步算

出

为了分析计算公式的精确度，通常可用泰勒展开将()1n y x +在n x 处展开，则有

()()()()()()2

''11,,.2

++'=+=++∈n n n n n n n n h y x y x h y x y x h y x x ξξ 在()n n y y x =的前提下，()()()(),,.n n n n n f x y f x y x y x '==可得欧拉法（3）的误差为 容易看出，欧拉法（3）式具有一阶精度.

2）向后欧拉方法：如果对微分方程（1）从n x 到1n x +积分，得

()()()()1

1,n n x n n x y x y x f t y t dx ++=+?，

（4） 如果（4）式右端积分用右矩形公式()()11,n n h f x y x ++?近似，则得到另一个公式 ()111,n n n n y y hf x y +++=+， （5） 称为后退欧拉法.

值得一提的是:后退欧拉法与欧拉公式有着本质的区别，后者是关于1n y +的直接计算公式，它是显式的，而（5）式的右端含有关于1n y +的表达式，它是隐式的.在利用后退欧拉法时，我们通常利用迭代法求解，实质就是逐步显示化.具体迭代过程如下：

首先利用欧拉公式(0)1(,)+=+?n n n n y y h f x y 给出迭代初值(0)1+n y ，把它代入（5）式的

右端，使之转化为显式，直接计算得11

(1)(0)1(,)+++=+?n n n n y y h f x y .如此反复进行，得 (1)()111(,)++++=+?k k n n n n y y h f x y 0,1,k =L ，

则得到后退欧拉法的迭代公式

(0)1(1)()111(,)(,)

+++++?=+??=+??n n n n k k n n n n y y h f x y y y h f x y ， 可以看出，后退欧拉法具有一阶精度，且计算比较麻烦.

1.2.2梯形方法

为得到比欧拉法精确度高的计算公式，在等式（4）式右端积分中若用梯形求积公式近似，并用n y 代替()n y x ，1n y +代替()1n y x +，则得

()()111,,2

n n n n n n h y y f x y f x y +++=++????， （6） 称其为梯形方法.

梯形方法与后退欧拉法一样，都是隐式单步法，可用迭代法求解，其迭代公式为 ()()(0)1(1)()111(,),,2

+++++?=+?????=++????n n n n k k n n n n n n y y h f x y h y y f x y f x y . （7） 为了分析梯形公式的收敛性，将（6）与（7）式相减，得

()()(1)()111111,,2

k k n n n n n n h y y f x y f x y +++++++??-=

-??，0,1,2,k =L 因为(),f x y 满足Lipschitz 条件，于是有(1)()11112

+++++-≤-k k n n n n hL y y y y ，其中L 为(),f x y 关于y 的Lipschitz 常数.如果选取h 充分小，使得12hL <，则当k →∞时有(1)11

+++→k n n y y ，这说明迭代过程（7）式是收敛的[4].容易推导得出梯形法（7）式是二阶方法.

经分析，梯形方法虽然提高了精度，但是以增加计算量为代价的.从上述的迭代公式可以看出，每迭代一次都要重新计算(),f x y 的值，而且迭代又要进行若干次，计算相当的复杂.为此，有没有比较简便的计算方法呢？下面给出改进的欧拉方法.

1.2.3改进的欧拉方法

由前面的讨论可知，梯形法计算相对复杂，现对上面的梯形法进行简化，具体方法是只计算一两次就转入下一步的计算，先用欧拉公式（3）求得一个初步的近似解

1n y +，称为预测值，再利用公式（6）把它校正一次，这样建立的预测-校正系统通常

称为改进的欧拉公式.具体公式如下

()()()11,,,2n n n n n n n n h y y f x y f x y hf x y ++??=+++?

? （8）改进的欧拉法与梯形法一样，是二阶方法.

1.2.4 Runge-Kutta 方法

由前面讨论可知，从（4）式

可以看出，若要使公式阶数提高，就必须使右端积分的数值求积公式精度提高，它必然要增加求积积点，为此将（4）式的右端用求积公式表示为

()()()()1

1,,n n r

x i n i n i x i f x y x dx h c f x h y x h λλ+=≈++∑?， （9） 一般来说，点数r 越多，精度越高，上式右端相当于增量函数(),,x y h ?，为得到便于计算的显式方法，将公式（9）表示为：

()1,,,n n n n y y h x y h ?+=+

（10） 其中

()()1111,,,,,2,r n n i i i n n i i n i n ij j j x y h c K K f x y K f x h y h K i r ?λμ=-=??=??=?????=++= ?????

∑∑L （11） 这里,,i i ij c λμ均为常数. i c 为加权因子，i K 为第i 段斜率，共有r 段.我们把（10）和

（11）称为r 级显式Runge-Kutta 法，简称为R-K 方法.下面给出其中最经典最常用的一个公式：

()()()11234121324322,6,,,22,,22,.n n n n n n n n n n h y y K K K K K f x y h h K f x y K h h K f x y K K f x h y hK +?=++++??=?????=++ ????????=++ ??????=++?

（12） Runge-Kutta 方法作为一种重要的单步方法，具有很高的实用价值，它关于初值是稳定的，其解连续地依赖于初值，是一类便于应用的单步法，为了计算1n y +，只用到前面一步的值n y 即可，因此每步的步长可以独立取定.常用的Runge-Kutta 方法精度较高，为了达到预定的精度，与欧拉方法与梯形法相比，步长h 可取得大些，求解区间上的总步数可以少些.但Runge-Kutta 方法也有些缺点，比如四阶Runge-Kutta 方法每算一步需要四次计算(),f x y 的值，计算量较大（对于复杂的

(),f x y 而言）

． 2 数值方法的应用实例[5-9]

例1 对于初值问题()1001

y y y '=-???=??，分别用欧拉法、改进的欧拉法，梯形法求()1y 的近似值.

解：易得该方程的解析解()10x y x e -=，()1 4.5400e-005y =，为比较，将按不同数值计算方法所得结果列表如下：

表 1 三种不同方法的数值结果

图 1 三种不同方法数值解与精确解的误差曲线从表1中可以看出：当0.2

h=时，三种方法均不稳定，计算结果严重偏离精确值；h=时，改进后的欧拉和梯形法均稳定，但欧拉法效果很差；当0.01

h≤时，三种

0.1

方法均稳定，但精确度有区别．可以看出，h越小，计算结果越好，要想计算结果充分接近于解析解还须取较小的h值．

图1反映的步长0.01

h=时，三种数值方法的所得数值解与解析解在[0,1]区间的误差曲线，由图可知，在步长相同的情况下，梯形法的精确度略高于改进的欧拉法；改进的欧拉法和梯形法精确度都明显高于欧拉法．

例2用欧拉法、改进的欧拉法和Runge-Kutta法求解初值问题

并比较三种方法的结果．

解：方程为1

n=-的伯努利方程，可求得解析解为

现用MATLAB软件编程，用题目要求的方法求解，可得如下图示结果：

图2 （a）步长为0.2时R-K法和解析解比较

图2 （b）步长为0.2时改进的Euler法和解析解比较

图2 （c）步长为0.2时欧拉法和解析解比较

上图2(a)，(b)，(c)描述的是步长为0.2时，用欧拉法、改进的欧拉法，Runge-Kutta 法求解方程所得的数值解与解析解之间的对比图．由图可知，Runge-Kutta法所得数值解曲线和解析解曲线吻合的很好，改进的欧拉法和欧拉法随着计算的进行，数值解和解析解之间误差逐步增大，但改进的欧拉法效果要好于欧拉法．

图3 (a) 步长为0.1时Euler法和解析解比较图3 (b) 步长为0.1时改进的Euler法和解析解比较

图3 (c) 步长为0.1时Runge-Kutta法和解析解比较

上图3 (a)，(b)，(c)描述的是步长为0.1时，用欧拉法、改进的欧拉法，Runge-Kutta

法求解方程所得的数值解与解析解之间的对比图．由图可知，改进的欧拉法和Runge-Kutta法所得数值解曲线和解析解曲线吻合的很好，而欧拉法随着计算的进行数值解和解析解之间误差逐步增大．

相应的程序如下：

主程序

x=0:0.2:1;

jxj=exp(2*x).*(1./exp(4*x) + (2*x)./exp(4*x)).^(1/2);

y=Euler(@ff,0,1,0.2,1);gy=geuler(@ff,0,1,0.2,1);Ry=RK(@ff,0,1,0.2,1); figure(1);plot(x,jxj,x,Ry,'*');figure(2);plot(x,jxj,x,gy,'*');figure(3);p lot(x,jxj,x,y,'*')

欧拉法程序

function y=Euler(f,a,b,h,y0)

n=(b-a)/h;x=a:h:b;y=zeros(n+1,1);y(1)=y0;

for i=1:n

y(i+1)=y(i)+h*feval(f,x(i),y(i));

end

改进的欧拉法程序

function gy=geuler(f,a,b,h,y0)

n=(b-a)/h;x=a:h:b;y=zeros(n+1,1);y(1)=y0;

for i=1:n

yp=y(i)+h*feval(f,x(i),y(i));

yc=y(i)+h*feval(f,x(i+1),yp);

y(i+1)=(yp+yc)/2;

end

gy=y;

Runge-Kutta法程序

function Ry=RK(f,a,b,h,y0)

n=(b-a)/h;x=a:h:b;y=zeros(n+1,1);y(1)=y0;

for i=1:n

k1=feval(f,x(i),y(i));

k2=feval(f,x(i)+h/2,y(i)+h*k1/2);

k3=feval(f,x(i)+h/2,y(i)+h*k2/2);

k4=feval(f,x(i+1),y(i)+h*k3);

y(i+1)=y(i)+h*(k1+2*k2+2*k3+k4)/6;

end

Ry=y;

3 微分方程数值解法在实际生活中的应用

3.1应用实例:耐用消费新产品的销售规律模型

一种新产品进入市场以后，常常会经历销售量首先慢慢增加然后逐渐慢慢下降的一个过程，由此给出的时间—销售坐标系下的曲线称为产品的生命曲线,它的形状是钟形的.不过对于较耐用的消费品，情况会有所不一样，它的生命曲线会在开始有个小小的高峰，之后是段比较平坦的曲线，先下降，之后再上升，然后达到高峰，因此它是双峰型的曲线.

如何理解这种和传统产品的生命曲线理论相冲突的现象？澳大利亚学者斯蒂芬斯与莫赛观察到的购买耐用消费产品的人大概可分为两类：其一是非常善于接受新的事物，称其为“创新型”消费者，他们会经常从产品广告，制造商给出产品的说明书与商店样品中了解了产品功能与性能之后，再决定其否购买；其二是消费者相对保守，称其为“模仿型”消费者，他们往往会根据大部分已经购买产品的消费者实际使用的经验而提供的信息来决定是否购买产品，下面的例子经过分析，建立相应的数学模型，对这种现象给出了科学解释.

3.1.1 模型的建立

将顾客获得信息大致可分成两类，其一称之为“搜集型”，源自于产告说明、广告，“创新型”顾客获得这类消息后就可作出其是否购买；第二类称之为“体验型”，即消费者使用之后会获得真实体验，常常以口头的形式散播，“模仿型”类顾客会在获得这种信息之后才能决定是否购买.

设K 是潜在用户的总数，1K 与2K 分别为 “创新型”与“模仿型”的人数.再设()N t 是时刻t 已经购买的商品顾客的人数，而()1N t 与()2N t 分别为“创新型”与

“模仿型”的顾客的人数,再设()1A t 是时刻t 中已获得“搜集型”信息的人数，由于该部分的信息可直接由外部获得，或者可从已获得该类信息的人群中获得，因此类似巴斯模型，从而建立如下方程：

()()()()()()11112112,00,,0dA t K A t A t A dt

αααα=-+=> ， 获得“搜集型”信息的“创新型”消费者决定其是否购买的行为，满足如下方程： 对于“模仿型”的顾客，可从已经购买该产品“创新型”或者“模仿型”的顾客中获取信息，于是有

在这里，忽略消费者购买产品后需一段短暂的使用后才会散播体验信息的滞后作用． 综上，斯蒂芬斯—莫赛模型是一常微分方程组的初值问题：

而()()()12N t N t N t =+为时刻t 购买该商品的总人数.

3.1.2 模型的求解

对于斯蒂芬斯—莫赛模型中()2N t 的解析解则不能求出，于是可以用四阶Runge Kutta -公式求得，且在它的精度要求达到很高情形下求出()2N t .

用MATLAB 软件求解，相应的程序及结果如下

function RK=RKFC(fc,a,b,h,y0)

n=(b-a)/h;x=a:h:b;m=length(y0);y=zeros(n+1,m);y(1,:)=y0;

for i=1:n

k1=feval(fc,x(i),y(i,:));k2=feval(fc,x(i)+h/2,y(i,:)+h*k1/2);

k3=feval(fc,x(i)+h/2,y(i,:)+h*k2/2);k4=feval(fc,x(i+1),y(i,:)+h*k3); y(i+1,:)=y(i,:)+h*(k1+2*k2+2*k3+k4)/6;

end

RK=y;

function f=FC(x,y)

k1=50; k2=70;al=0.0013;be=0.0013;ga=0.0015;

f=[(k1-y(1))*(al+be*y(1)),ga*(k2-y(2))*(y(1)+y(2))];

x=0:0.3:24;

RK=RKFC(@FC,0,24,0.3,[0,0])

figure(1)；plot(x,RK(:,1),'+',x,RK(:,2),'*')；legend('N1(t)','N2(t)',2) figure(2)；plot(x,RK(:,1)+RK(:,2),'+')；legend('N(t)',2)

图 4 ()()12,N t N t 与时间关系图

图 5 ()N t 与[0,25]时间段关系图 由此例可以看出，微分方程数值解在实际生活有着广泛的应用，而数值解法中的Runge-Kutta 方法是一种很重要且应用很广泛的算法.

结语

微分方程数值解是求解微分方程的一种很重要且应用范围很广的方法，而常用的数值解法如欧拉法、改进的欧拉法、梯形法和Runge-Kutta 方法在一定程度上都有自己的优缺点，理解和掌握各种方法的应用范围，用MATLAB 编制各种方法求解实际问题的通用程序，对用微分方程数值解理论解决现实生活中的实际问题有很重要的意义．

参考文献

[1] 李庆扬,王能超,易大义.数值分析[M].北京:清华大学出版社,清华大学出版社,2001.

[2] 冯康.数值计算方法[M].杭州:浙江大学出版社,2003.

[3] 封建湖.计算方法典型题分析解集[M].西安:西北工业大学出版社,2000.

[4] 胡建伟,汤怀民.微分方程数值解法(第二版)[M].北京:科学出版社,2007.

[5] 王能超.计算方法简明教程[M].北京:高等教育出版社,2004.

[6] Nash S G.A history of scientific computing.[M] New York:ACM Press,1990.

[7] 于丽妮. ODE 问题解析解及数值解的matlab 实现[J]. 电脑知识与技术.2012,8(14):3303-3305.

[8] 霍晓成. 常微分方程数值解法的研究[J]. 临沂师范学院学报. 2011,(6):19-23.

[6] 王国立, 陈 瑛.非线性微分方程迭代算法及其在物理学中的应用[J]. 长春师范学院学报( 自然科学版). 2006, 25(2):10-12.