微分方程数值解讲解Word文档下载推荐.docx
《微分方程数值解讲解Word文档下载推荐.docx》由会员分享，可在线阅读，更多相关《微分方程数值解讲解Word文档下载推荐.docx（15页珍藏版）》请在冰豆网上搜索。
其中多元函数与fx,y有关，当含有yn1时，方法是隐式的；
若中不含yn1，则为显式方法，所以显式单步法可表示为
yn1ynh（xn,yn,h）.
（2）
设yx是初值问题
（1）的准确解，称Tn1yxn1yxnh（xn,yxn,h）为
显式单步法
（2）的局部截断误差.若存在最大正整数p,使显式单步法
（2）式的局
部截断误差满足Tn1yxhyxhx,y,hOhp1，则称
（2）式有p阶精度.
1.2几种常用的数值解法及其分析、比较
1.2.1欧拉法与后退欧拉法
1）欧拉法：
欧拉曾简单地用差分代替微分，即利用公式
将初值问题
（1）离散化，则问题
（1）可化为
yn1ynhf（xn,yn）,xnx0nh，（3）
此方法称为欧拉法.
欧拉方法的几何意义在数值计算公式中体现了出来.在xy平面上，一阶微分方程的解yyx称作它的积分曲线.积分曲线上一点x,y的切线斜率等于函数fx,y，按函数fx,y在xy平面上建立一个方向场，那么，积分曲线上每一点的切线方向均与方向场在该点的方向相一致.
基于上述几何解释，从初始点P0（x0,y0）出发，先依方向场在该点的方向上推进到xx1上一点P1，再从P1依方向场的方向推进到xx2上一点P2，循环前进便作出一条折线P0P1P2，因此欧拉方法又称为折线法.若初值y0已知，则由（3）式可逐步算出
y1y0hf（x0,y0）,
y2y1hf（x1,y1）,
为了分析计算公式的精确度，通常可用泰勒展开将yxn1在xn处展开，则有
nxn,xn1.
yxn1yxnhyxnyxnhh2y'
'
n,
在ynyxn的前提下，fxn,ynf
xn,yxnyxn.可得欧拉法（3）的误差为
h2h2
yxn1yn1h2ynh2yxn.
容易看出，欧拉法（3）式具有一阶精度.
2）向后欧拉方法：
如果对微分方程
（1）从xn到xn1积分，得
4）
xn1
yxn1yxnxft,ytdx，
xn
如果（4）式右端积分用右矩形公式hfxn1,yxn1近似，则得到另一个公式
称为后退欧拉法.
值得一提的是:
后退欧拉法与欧拉公式有着本质的区别，后者是关于yn1的直接
计算公式，它是显式的，而（5）式的右端含有关于yn1的表达式，它是隐式的.在利用后退欧拉法时，我们通常利用迭代法求解，实质就是逐步显示化.具体迭代过程如下：
首先利用欧拉公式yn（0）1ynhf（xn,yn）给出迭代初值yn（0）1，把它代入（5）式的
右端，使之转化为显式，直接计算得y
（1）ynhf（xn1,y（0））.如此反复进行，得
y（nk11）ynhf（xn1,y（nk）1）k0,1,，
则得到后退欧拉法的迭代公式
yn（0）1ynhf（xn,yn），
yn（k11）ynhf（xn1,yn（k）1）
可以看出，后退欧拉法具有一阶精度，且计算比较麻烦.
1.2.2梯形方法
为得到比欧拉法精确度高的计算公式，在等式（4）式右端积分中若用梯形求积
公式近似，并用yn代替yxn，yn1代替yxn1，则得
6）
yn1ynh2fxn,ynfxn1,yn1，称其为梯形方法.
梯形方法与后退欧拉法一样，都是隐式单步法，可用迭代法求解，其迭代公式为
7）
yn（0）1ynhf（xn,yn）
yn（k11）ynh2fxn,ynfxn1,yn（k）1
为了分析梯形公式的收敛性，将（6）与（7）式相减，得
yn1yn（k11）h2fxn1,yn1fxn1,yn（k1），k0,1,2,
关于y的Lipschitz常数.如果选取h充分小，使得hL1，则当k时有
2
y（nk11）yn1，这说明迭代过程（7）式是收敛的[4].容易推导得出梯形法（7）式是二
阶方法.
经分析，梯形方法虽然提高了精度，但是以增加计算量为代价的.从上述的迭代
公式可以看出，每迭代一次都要重新计算fx,y的值，而且迭代又要进行若干次，计算相当的复杂.为此，有没有比较简便的计算方法呢？
下面给出改进的欧拉方法.
1.2.3改进的欧拉方法
由前面的讨论可知，梯形法计算相对复杂，现对上面的梯形法进行简化，具体方法是只计算一两次就转入下一步的计算，先用欧拉公式（3）求得一个初步的近似解yn1，称为预测值，再利用公式（6）把它校正一次，这样建立的预测-校正系统通常称为改进的欧拉公式.具体公式如下
h
yn1ynfxn,ynfxn1,ynhfxn,yn（8）
改进的欧拉法与梯形法一样，是二阶方法.
1.2.4Runge-Kutta方法
由前面讨论可知，从（4）式
yxn1yxnxfx,yxdx
可以看出，若要使公式阶数提高，就必须使右端积分的数值求积公式精度提高，它必然要增加求积积点，为此将（4）式的右端用求积公式表示为
xn1r
xfx,yxdxhcifxnih,yxnih，（9）
xni1inini
一般来说，点数r越多，精度越高，上式右端相当于增量函数x,y,h，为得到便于计算的显式方法，将公式（9）表示为：
yn1ynhxn,yn,h,（10）
其中
r
xn,yn,hciKi
i1
K1fxn,yn（11）
Kifxnih,ynhijKj,i2,r
j1
这里ci,i,ij均为常数.ci为加权因子，Ki为第i段斜率，共有r段.我们把（10）和
（11）称为r级显式Runge-Kutta法，简称为R-K方法.下面给出其中最经典最常用的一个公式：
yn1ynK12K22K3K4,
6
hh
K3fxnh2,ynh2K2
K4fxnh,ynhK3.
Runge-Kutta方法作为一种重要的单步方法，具有很高的实用价值，它关于初值
是稳定的，其解连续地依赖于初值，是一类便于应用的单步法，为了计算yn1，只用到前面一步的值yn即可，因此每步的步长可以独立取定.常用的Runge-Kutta方法精度较高，为了达到预定的精度，与欧拉方法与梯形法相比，步长h可取得大些，求解区间上的总步数可以少些.但Runge-Kutta方法也有些缺点，比如四阶Runge-Kutta方法每算一步需要四次计算fx,y的值，计算量较大（对于复杂的fx,y而言）．
2数值方法的应用实例[5-9]
y10y
例1对于初值问题，分别用欧拉法、改进的欧拉法，梯形法求y1的y01
近似值.
解：
易得该方程的解析解yxe10x，y14.5400e-005，为比较，将按不同数值计算方法所得结果列表如下：
欧拉法
改进的欧拉法
梯形法
0.2
-1
1
0.1
9.7656E-004
5.4994E-005
0.01
2.6561E-005
4.6223E-005
4.5026E-005
0.001
4.3717E-005
4.5408E-005
4.5396E-005
0.0001
4.5173E-005
4.5400E-005
4.5400E-005
表1三种不同方法的数值结果
0.02
0.015
0.005
欧拉法误差改进欧拉法误差梯形法误差
图1三种不同方法数值解与精确解的误差曲线
从表1中可以看出：
当h0.2时，三种方法均不稳定，计算结果严重偏离精确值；
h0.1时，改进后的欧拉和梯形法均稳定，但欧拉法效果很差；
当h0.01时，三种方法均稳定，但精确度有区别．可以看出，h越小，计算结果越好，要想计算结果充分接近于解析解还须取较小的h值．
图1反映的步长h0.01时，三种数值方法的所得数值解与解析解在[0,1]区间的误差曲线，由图可知，在步长相同的情况下，梯形法的精确度略高于改进的欧拉法；
改进的欧拉法和梯形法精确度都明显高于欧拉法．
例2用欧拉法、改进的欧拉法和Runge-Kutta法求解初值问题
dyy2x,x0,1
dxy
y01
并比较三种方法的结果．
方程为n1的伯努利方程，可求得解析解为
ye2xe4x2xe4x
现用MATLAB软件编程，用题目要求的方法求解，可得如下图示结果：
图2（a）步长为0.2时R-K法和解析解比较
1.8
1.6
1.4
1.2
解析解
改进的Euler法
图2（b）步长为0.2时改进的Euler法和解析解比较
图2（c）步长为0.2时欧拉法和解析解比较
上图2（a），（b），（c）描述的是步长为0.2时，用欧拉法、改进的欧拉法，Runge-Kutta法求解方程所得的数值解与解析解之间的对比图．由图可知，Runge-Kutta法所得数值解曲线和解析解曲线吻合的很好，改进的欧拉法和欧拉法随着计算的进行，数值解
和解析解之间误差逐步增大，但改进的欧拉法效果要好于欧拉法．
解析解Euler法
Euler法
00.20.40.60.81
图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）;
plot（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:
n
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）;
10k4=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已经购买的商品顾客的人数，而N1t与N2t分别为“创新型”与“模仿型”的顾客的人数,再设A1t是时刻t中已获得“搜集型”信息的人数，由于该部分的信息可直接由外部获得，或者可从已获得该类信息的人群中获得，因此类似巴斯模型，从而建立如下方程：
11
dA1tK1A1t12A1t,A00,1,20，
dt
获得“搜集型”信息的“创新型”消费者决定其是否购买的行为，满足如下方程：
dN1tK1N1tN1t,N100,,0
对于“模仿型”的顾客，可从已经购买该产品“创新型”或者“模仿型”的顾客中获取信息，于是有
在这里，忽略消费者购买产品后需一段短暂的使用后才会散播体验信息的滞后作用．
综上，斯蒂芬斯—莫赛模型是一常微分方程组的初值问题：
dN2t
2K2N2tN1tN2t,dt
N100,N2（0）0
而NtN1tN2t为时刻t购买该商品的总人数.
3.1.2模型的求解
对于斯蒂芬斯—莫赛模型中N2t的解析解则不能求出，于是可以用四阶
RungeKutta公式求得，且在它的精度要求达到很高情形下求出N2t.
用MATLAB软件求解，相应的程序及结果如下
functionRK=RKFC（fc,a,b,h,y0）n=（b-a）/h;
m=length（y0）;
y=zeros（n+1,m）;
y（1,:
）=y0;
nk1=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;
12
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）；
plot（x,RK（:
1）+RK（:
N（t）'
2）
8
图4N1t,N2t与时间关系图
图5Nt与[0,25]时间段关系图
由此例可以看出，微分方程数值解在实际生活有着广泛的应用，而数值解法中的Runge-Kutta方法是一种很重要且应用很广泛的算法.
13
结语
微分方程数值解是求解微分方程的一种很重要且应用范围很广的方法，而常用的数值解法如欧拉法、改进的欧拉法、梯形法和Runge-Kutta方法在一定程度上都有自己的优缺点，理解和掌握各种方法的应用范围，用MATLAB编制各种方法求解实际问题的通用程序，对用微分方程数值解理论解决现实生活中的实际问题有很重要的意义．
14
参考文献
[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-3305.
[8]霍晓成.常微分方程数值解法的研究[J].临沂师范学院学报.2011,（6）:
19-23.
自然科学版）.
[6]王国立,陈瑛.非线性微分方程迭代算法及其在物理学中的应用[J].长春师范学院学报（
2006,25
（2）:
10-12.
15
升级会员 