微分方程数值解docxWord文档下载推荐.docx
《微分方程数值解docxWord文档下载推荐.docx》由会员分享，可在线阅读，更多相关《微分方程数值解docxWord文档下载推荐.docx（21页珍藏版）》请在冰豆网上搜索。
解问题．求解微分方程虽然有多种解析方法，但根据工程和科学实践问题所得到的微
分方程往往很复杂，在很多情况下不能或很难给出解析解，有时即使能求出形式解，
也往往因形式过于复杂或计算量太大而不实用，因此从实际问题中归结出来的微分方
程主要依靠数值解法．
定义2微分方程数值解：
对初值问题
（1）寻求数值解就是寻求解yx在一
系列离散节点
上的近似解
y0,y1,y2,L,yn,yn1,L，相邻两个节点的间距hnxn1xn称为步长.在一
般情况下假定hihi0,1,L为常数，这时节点为xnx0nh,n0,1,2,L．
要求微分方程数值解，首先要建立数值算法，即对初值问题
（1）中的方程离散
化，建立求解数值解法的递推公式．一类是计算yn1时只用到前一点的值yn，称为单
步法；
另一类是用到
y前面k点的值
n1
yyLy称为k步法.
,,
nn1nk1
对初值问题
（1）式的单步法可用一般形式表示为
y1yh（x,y,y1,h），
nnnnn
其中多元函数与fx,y有关，当含有
y时，方法是隐式的；
若中不含yn1，
则为显式方法，所以显式单步法可表示为
yyhxyh.
nnnn
1（,,）
（2）
设yx是初值问题
（1）的准确解，称
Tyxyxhxyxh为
11（,,）
显式单步法
（2）的局部截断误差.若存在最大正整数p,使显式单步法
（2）式的局
部截断误差满足
p1
T1yxhyxhx,y,hOh，则称
（2）式有p阶精度.
n
1.2几种常用的数值解法及其分析、比较
1.2.1欧拉法与后退欧拉法
1）欧拉法：
欧拉曾简单地用差分代替微分，即利用公式
将初值问题
（1）离散化，则问题
（1）可化为
y1yhf（x,y）,xnx0nh，（3）
此方法称为欧拉法.
欧拉方法的几何意义在数值计算公式中体现了出来.在xy平面上，一阶微分方
程的解yyx称作它的积分曲线.积分曲线上一点x,y的切线斜率等于函数
fx,y，按函数fx,y在xy平面上建立一个方向场，那么，积分曲线上每一点的切
线方向均与方向场在该点的方向相一致.
基于上述几何解释，从初始点P0（x0,y0）出发，先依方向场在该点的方向上推进到
xx上一点P1，再从P1依方向场的方向推进到xx2上一点P2，循环前进便作出一
1
条折线
PPPL，因此欧拉方法又称为折线法.若初值y0已知，则由（3）式可逐步算
012
出
为了分析计算公式的精确度，通常可用泰勒展开将
yx在xn处展开，则有
2
h
'
yx1yxhyxyxhy,x,x1.
nnnnnnnn
在ynyxn的前提下，fxn,ynfxn,yxnyxn.可得欧拉法（3）的误差为
容易看出，欧拉法（3）式具有一阶精度.
2）向后欧拉方法：
如果对微分方程
（1）从xn到xn1积分，得
x
yx1yxft,ytdx，（4）
nn
如果（4）式右端积分用右矩形公式
hfx1,yx1近似，则得到另一个公式
yyhfxy，（5）
11,1nnnn
称为后退欧拉法.
值得一提的是:
后退欧拉法与欧拉公式有着本质的区别，后者是关于yn1的直接
计算公式，它是显式的，而（5）式的右端含有关于
y的表达式，它是隐式的.在利
用后退欧拉法时，我们通常利用迭代法求解，实质就是逐步显示化.具体迭代过程如
下：
首先利用欧拉公式
（0）
y1yhf（x,y）给出迭代初值
y，把它代入（5）式的
右端，使之转化为显式，直接计算得
（1）（0）
yyhf（x,y）.如此反复进行，得
nnn1n
11
（k1）（k）
y1yhf（x1,y1）k0,1,L，
则得到后退欧拉法的迭代公式
yyhf（x,y）
n1nnn
（k1）（k）
n1nn1n1
，
可以看出，后退欧拉法具有一阶精度，且计算比较麻烦.
1.2.2梯形方法
为得到比欧拉法精确度高的计算公式，在等式（4）式右端积分中若用梯形求积
公式近似，并用
y代替
yx，yn1代替
yx，则得
yyfxyfxy，（6）
nnnnnn
1,1,1
称其为梯形方法.
梯形方法与后退欧拉法一样，都是隐式单步法，可用迭代法求解，其迭代公式为
（k1）（k）
yyfx,yfx,y
n1nnnn1n1
.（7）
为了分析梯形公式的收敛性，将（6）与（7）式相减，得
（k1）（k）
yyfxyfxy，k0,1,2,L
111,11,1nnnnnn
因为fx,y满足Lipschitz条件，于是有
hL
yyyy，其中L为fx,y
11211
关于y的Lipschitz常数.如果选取h充分小，使得1
，则当k时有
（k1）
yy，这说明迭代过程（7）式是收敛的
n1n1
[4].容易推导得出梯形法（7）式是二
阶方法.
经分析，梯形方法虽然提高了精度，但是以增加计算量为代价的.从上述的迭代
公式可以看出，每迭代一次都要重新计算fx,y的值，而且迭代又要进行若干次，
计算相当的复杂.为此，有没有比较简便的计算方法呢？
下面给出改进的欧拉方法.
1.2.3改进的欧拉方法
由前面的讨论可知，梯形法计算相对复杂，现对上面的梯形法进行简化，具体方
法是只计算一两次就转入下一步的计算，先用欧拉公式（3）求得一个初步的近似解
y，称为预测值，再利用公式（6）把它校正一次，这样建立的预测-校正系统通常
称为改进的欧拉公式.具体公式如下
yyfxyfxyhfxy（8）
1,1,,
改进的欧拉法与梯形法一样，是二阶方法.
1.2.4Runge-Kutta方法
由前面讨论可知，从（4）式
可以看出，若要使公式阶数提高，就必须使右端积分的数值求积公式精度提高，它必
然要增加求积积点，为此将（4）式的右端用求积公式表示为
r
fxyxdxhcfxhyxh，（9）
,
inini
i
一般来说，点数r越多，精度越高，上式右端相当于增量函数x,y,h，为得到便于
计算的显式方法，将公式（9）表示为：
yyhxyh（10）
1,,,nnnn
其中
x,y,hcK
nnii
i1
Kfx,y
1nn
（11）
Kfxh,yhK,i2,Lr
ininijj
j1
这里ci,i,ij均为常数.
c为加权因子，Ki为第i段斜率，共有r段.我们把（10）和
（11）称为r级显式Runge-Kutta法，简称为R-K方法.下面给出其中最经典最常用
的一个公式：
yyK2K2KK,
n1n1234
6
Kfx,y,
hh
Kfx,yK
2nn1
22
（12）
Kfx,yK,
3nn2
Kfxh,yhK.
4nn3
Runge-Kutta方法作为一种重要的单步方法，具有很高的实用价值，它关于初值
是稳定的，其解连续地依赖于初值，是一类便于应用的单步法，为了计算yn1，只用
到前面一步的值
y即可，因此每步的步长可以独立取定.常用的Runge-Kutta方
法精度较高，为了达到预定的精度，与欧拉方法与梯形法相比，步长h可取得大些，
求解区间上的总步数可以少些.但Runge-Kutta方法也有些缺点，比如四阶
Runge-Kutta方法每算一步需要四次计算fx,y的值，计算量较大（对于复杂的
fx,y而言）．
[5-9]2数值方法的应用实例
例1对于初值问题
y10y
y01
，分别用欧拉法、改进的欧拉法，梯形法求y1的
近似值.
解：
易得该方程的解析解
10x
yxe，y14.5400e-005，为比较，将按不同
数值计算方法所得结果列表如下：
表1三种不同方法的数值结果
欧拉法改进的欧拉法梯形法
0.2-111
0.109.7656E-0045.4994E-005
0.012.6561E-0054.6223E-0054.5026E-005
0.0014.3717E-0054.5408E-0054.5396E-005
0.00014.5173E-0054.5400E-0054.5400E-005
图1三种不同方法数值解与精确解的误差曲线
从表1中可以看出：
当h0.2时，三种方法均不稳定，计算结果严重偏离精确值；
h0.1时，改进后的欧拉和梯形法均稳定，但欧拉法效果很差；
当h0.01时，三种
方法均稳定，但精确度有区别．可以看出，h越小，计算结果越好，要想计算结果充
分接近于解析解还须取较小的h值．
图1反映的步长h0.01时，三种数值方法的所得数值解与解析解在[0,1]区间的误差
曲线，由图可知，在步长相同的情况下，梯形法的精确度略高于改进的欧拉法；
改进
的欧拉法和梯形法精确度都明显高于欧拉法．
例2用欧拉法、改进的欧拉法和Runge-Kutta法求解初值问题
并比较三种方法的结果．
方程为n1的伯努利方程，可求得解析解为
现用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,'
）
欧拉法程序
functiony=Euler（f,a,b,h,y0）
n=（b-a）/h;
x=a:
h:
b;
y=zeros（n+1,1）;
y
（1）=y0;
fori=1:
y（i+1）=y（i）+h*feval（f,x（i）,y（i））;
end
改进的欧拉法程序
functiongy=geuler（f,a,b,h,y0）
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;
gy=y;
Runge-Kutta法程序
functionRy=RK（f,a,b,h,y0）
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;
Ry=y;
3微分方程数值解法在实际生活中的应用
3.1应用实例:
耐用消费新产品的销售规律模型
一种新产品进入市场以后，常常会经历销售量首先慢慢增加然后逐渐慢慢下降的
一个过程，由此给出的时间—销售坐标系下的曲线称为产品的生命曲线,它的形状是
钟形的.不过对于较耐用的消费品，情况会有所不一样，它的生命曲线会在开始有个
小小的高峰，之后是段比较平坦的曲线，先下降，之后再上升，然后达到高峰，因此
它是双峰型的曲线.
如何理解这种和传统产品的生命曲线理论相冲突的现象？
澳大利亚学者斯蒂芬
斯与莫赛观察到的购买耐用消费产品的人大概可分为两类：
其一是非常善于接受新的
事物，称其为“创新型”消费者，他们会经常从产品广告，制造商给出产品的说明书
与商店样品中了解了产品功能与性能之后，再决定其否购买；
其二是消费者相对保守，
称其为“模仿型”消费者，他们往往会根据大部分已经购买产品的消费者实际使用的
经验而提供的信息来决定是否购买产品，下面的例子经过分析，建立相应的数学模型，
对这种现象给出了科学解释.
3.1.1模型的建立
将顾客获得信息大致可分成两类，其一称之为“搜集型”，源自于产告说明、广
告，“创新型”顾客获得这类消息后就可作出其是否购买；
第二类称之为“体验型”，
即消费者使用之后会获得真实体验，常常以口头的形式散播，“模仿型”类顾客会在
获得这种信息之后才能决定是否购买.
设K是潜在用户的总数，K1与K2分别为“创新型”与“模仿型”的人数.再
设Nt是时刻t已经购买的商品顾客的人数，而
Nt与N2t分别为“创新型”与
“模仿型”的顾客的人数,再设
At是时刻t中已获得“搜集型”信息的人数，由于
该部分的信息可直接由外部获得，或者可从已获得该类信息的人群中获得，因此类似
巴斯模型，从而建立如下方程：
dAt
dt
K1A1t12A1t,A00,1,20，
获得“搜集型”信息的“创新型”消费者决定其是否购买的行为，满足如下方程：
对于“模仿型”的顾客，可从已经购买该产品“创新型”或者“模仿型”的顾客中获
取信息，于是有
在这里，忽略消费者购买产品后需一段短暂的使用后才会散播体验信息的滞后作用．
综上，斯蒂芬斯—莫赛模型是一常微分方程组的初值问题：
而
NtNtNt为时刻t购买该商品的总人数.
3.1.2模型的求解
对于斯蒂芬斯—莫赛模型中
Nt的解析解则不能求出，于是可以用四阶
RungeKutta公式求得，且在它的精度要求达到很高情形下求出N2t.
用MATLAB软件求解，相应的程序及结果如下
functionRK=RKFC（fc,a,b,h,y0）
m=length（y0）;
y=zeros（n+1,m）;
y（1,:
）=y0;
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;
RK=y;
functionf=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））];
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）；
1）+RK（:
N（t）'
图4
NtNt与时间关系图
1,2
图5Nt与[0,25]时间段关系图
由此例可以看出，微分方程数值解在实际生活有着广泛的应用，而数值解法中的
Runge-Kutta方法是一种很重要且应用很广泛的算法.
结语
微分方程数值解是求解微分方程的一种很重要且应用范围很广的方法，而常用的
数值解法如欧拉法、改进的欧拉法、梯形法和Runge-Kutta方法在一定程度上都有自
己的优缺点，理解和掌握各种方法的应用范围，用MATLAB编制各种方法求解实际问
题的通用程序，对用微分方程数值解理论解决现实生活中的实际问题有很重要的意
义．
参考文献
[1]李庆扬,王能超,易大义.数值分析[M].北京:
清华大学出版社,清华大学出版社,2001.
[2]冯康.数值计算方法[M].杭州:
浙江大学出版社,2003.
[3]封建湖.计算方法典型题分析解集[M].西安:
西北工业大学出版社,2000.
[4]胡建伟,汤怀民.微分方程数值解法（第二版）[M].北京:
科学出版社,2007.
[5]王能超.计算方法简明教程[M].北京:
高等教育出版社,2004.
[6]NashSG.Ahistoryofscientificcomputing.[M]NewYork:
ACMPress,1990.
[7]于丽妮.ODE问题解析解及数值解的matlab实现[J].电脑知识与技术.2012,8（14）:
3303-33
升级会员 