微分方程数值解docxWord文档下载推荐.docx

1、解问题求解微分方程虽然有多种解析方法， 但根据工程和科学实践问题所得到的微分方程往往很复杂，在很多情况下不能或很难给出解析解，有时即使能求出形式解，也往往因形式过于复杂或计算量太大而不实用， 因此从实际问题中归结出来的微分方程主要依靠数值解法定义 2 微分方程数值解：对初值问题（ 1）寻求数值解就是寻求解 y x 在一系列离散节点上的近似解y0 , y1, y2,L , yn , yn 1,L ，相邻两个节点的间距 hn xn 1 xn 称为步长 . 在一般情况下假定 hi h i 0,1,L 为常数，这时节点为 xn x0 nh, n 0,1, 2,L 要求微分方程数值解，首先要建立数值算法

2、，即对初值问题（ 1）中的方程离散化，建立求解数值解法的递推公式一类是计算 yn 1 时只用到前一点的值 yn ，称为单步法；另一类是用到y 前面k 点的值n 1y y L y 称为k 步法., , ,n n 1 n k 1对初值问题（ 1）式的单步法可用一般形式表示为y 1 y h （x , y , y 1, h） ，n n n n n其中多元函数 与 f x, y 有关，当 含有y 时，方法是隐式的；若 中不含 yn 1，则为显式方法，所以显式单步法可表示为y y h x y h .n n n n1 （ , , ）（2）设 y x 是初值问题（ 1）的准确解，称T y x y x h x

3、y x h 为1 1 （ , , ）显式单步法（ 2）的局部截断误差 . 若存在最大正整数 p , 使显式单步法（ 2）式的局部截断误差满足p 1T 1 y x h y x h x, y,h O h ，则称（2）式有 p 阶精度.n1.2 几种常用的数值解法及其分析、比较1.2.1 欧拉法与后退欧拉法1）欧拉法：欧拉曾简单地用差分代替微分，即利用公式将初值问题（ 1）离散化，则问题（ 1）可化为y 1 y h f （ x , y ）, xn x0 n h ， （3）此方法称为欧拉法 .欧拉方法的几何意义在数值计算公式中体现了出来 . 在xy 平面上，一阶微分方程的解 y y x 称作它的积分曲

4、线 . 积分曲线上一点 x, y 的切线斜率等于函数f x, y ，按函数 f x, y 在xy平面上建立一个方向场，那么，积分曲线上每一点的切线方向均与方向场在该点的方向相一致 .基于上述几何解释， 从初始点 P0（ x0, y0 ） 出发，先依方向场在该点的方向上推进到x x 上一点 P1，再从 P1 依方向场的方向推进到 x x2 上一点 P2 ，循环前进便作出一1条折线P PP L ，因此欧拉方法又称为折线法 . 若初值 y0 已知，则由（3）式可逐步算0 1 2出为了分析计算公式的精确度，通常可用泰勒展开将y x 在 xn 处展开，则有2hy x 1 y x h y x y x h

5、y , x , x 1 .n n n n n n n n在 yn y xn 的前提下， f xn , yn f xn, y xn y xn .可得欧拉法（ 3）的误差为容易看出，欧拉法（ 3）式具有一阶精度 .2）向后欧拉方法：如果对微分方程（ 1）从 xn 到xn 1积分，得xy x 1 y x f t, y t dx ， （4）n n如果（4）式右端积分用右矩形公式h f x 1, y x 1 近似，则得到另一个公式y y hf x y ， （5）1 1 , 1 n n n n称为后退欧拉法 .值得一提的是 : 后退欧拉法与欧拉公式有着本质的区别，后者是关于 yn 1 的直接计算公式，它是

6、显式的，而（ 5）式的右端含有关于y 的表达式，它是隐式的 . 在利用后退欧拉法时，我们通常利用迭代法求解，实质就是逐步显示化 . 具体迭代过程如下：首先利用欧拉公式（0）y 1 y h f （x , y ） 给出迭代初值y ，把它代入（ 5）式的右端，使之转化为显式，直接计算得（1） （0）y y h f （x , y ） . 如此反复进行，得n n n 1 n1 1（ k 1） （ k ）y 1 y h f （x 1, y 1） k 0,1,L ，则得到后退欧拉法的迭代公式y y h f （x , y ）n 1 n n n（ k 1） （k）n 1 n n 1 n 1，可以看出，后退欧拉法

7、具有一阶精度，且计算比较麻烦 .1.2.2 梯形方法为得到比欧拉法精确度高的计算公式，在等式（ 4）式右端积分中若用梯形求积公式近似，并用y 代替y x ，yn 1代替y x ，则得y y f x y f x y ， （6）n n n n n n1 , 1, 1称其为梯形方法 .梯形方法与后退欧拉法一样， 都是隐式单步法， 可用迭代法求解， 其迭代公式为（k 1） （ k）y y f x , y f x , yn 1 n n n n 1 n 1. （7）为了分析梯形公式的收敛性，将（ 6）与（7）式相减，得（k 1） （k ）y y f x y f x y ，k 0,1,2, L1 1 1,

8、1 1, 1 n n n n n n因为 f x, y 满足 Lipschitz 条件，于是有hLy y y y ，其中 L 为 f x,y1 1 2 1 1关于 y 的 Lipschitz 常数. 如 果选取 h 充分小，使得 1，则当 k 时有（ k 1）y y ，这说明迭代过程（ 7）式是收敛的n 1 n 14 . 容易推导得出梯形法（ 7）式是二阶方法.经分析，梯形方法虽然提高了精度，但是以增加计算量为代价的 . 从上述的迭代公式可以看出，每迭代一次都要重新计算 f x, y 的值，而且迭代又要进行若干次，计算相当的复杂 . 为此，有没有比较简便的计算方法呢？下面给出改进的欧拉方法 .

9、1.2.3 改进的欧拉方法由前面的讨论可知， 梯形法计算相对复杂， 现对上面的梯形法进行简化， 具体方法是只计算一两次就转入下一步的计算，先用欧拉公式（ 3）求得一个初步的近似解y ，称为预测值，再利用公式（ 6）把它校正一次，这样建立的预测 - 校正系统通常称为改进的欧拉公式 . 具体公式如下y y f x y f x y hf x y （8）1 , 1, ,改进的欧拉法与梯形法一样，是二阶方法 .1.2.4 Runge-Kutta 方法由前面讨论可知，从（ 4）式可以看出， 若要使公式阶数提高， 就必须使右端积分的数值求积公式精度提高， 它必然要增加求积积点，为此将（ 4）式的右端用求积公

10、式表示为rf x y x dx h c f x h y x h ， （9）, ,i n i n ii一般来说，点数 r 越多，精度越高，上式右端相当于增量函数 x,y,h ，为得到便于计算的显式方法，将公式（ 9）表示为：y y h x y h （10）1 , , , n n n n其中x , y ,h c Kn n i ii 1K f x ,y1 n n（11）K f x h, y h K ,i 2,L ri n i n ij jj 1这里 ci , i , ij 均为常数.c 为加权因子， Ki 为第 i 段斜率，共有 r 段. 我们把（10）和（11）称为 r 级显式 Runge-Kut

11、ta 法，简称为 R-K 方法. 下面给出其中最经典最常用的一个公式：y y K 2K 2K K ,n 1 n 1 2 3 46K f x , y ,h hK f x , y K2 n n 12 2（12）K f x , y K ,3 n n 2K f x h, y hK .4 n n 3Runge-Kutta 方法作为一种重要的单步方法，具有很高的实用价值，它关于初值是稳定的，其解连续地依赖于初值，是一类便于应用的单步法，为了计算 yn 1 ，只用到前面一步的值y 即可，因此每步的步长可以独立取定 . 常用的 Runge-Kutta 方法精度较高，为了达到预定的精度，与欧拉方法与梯形法相比，

12、步长 h可取得大些，求解 区间上的总步数可以少 些. 但 Runge-Kutta 方法也有些 缺点 ，比如 四阶Runge-Kutta 方法每算一步需要四次计算 f x, y 的值，计算量较大（对于复杂的f x, y 而言）5-9 2 数值方法的应用实例例 1 对于初值问题y 10 yy 0 1，分别用欧拉法、 改进的欧拉法， 梯形法求 y 1 的近似值.解：易得该方程的解析解10xy x e ， y 1 4.5400e-005，为比较，将按不同数值计算方法所得结果列表如下：表 1 三种不同方法的数值结果欧拉法 改进的欧拉法 梯形法0.2 -1 1 10.1 0 9.7656E-004 5.4

13、994E-0050.01 2.6561 E-005 4.6223 E -005 4.5026 E -0050.001 4.3717 E-005 4.5408 E -005 4.5396 E -0050.0001 4.5173 E-005 4.5400 E -005 4.5400 E -005图 1 三种不同方法数值解与精确解的误差曲线从表 1 中可以看出： 当h 0.2时，三种方法均不稳定， 计算结果严重偏离精确值；h 0.1时，改进后的欧拉和梯形法均稳定，但欧拉法效果很差；当 h 0.01时，三种方法均稳定，但精确度有区别可以看出， h 越小，计算结果越好，要想计算结果充分接近于解析解还须取

14、较小的 h值图 1 反映的步长 h 0.01时，三种数值方法的所得数值解与解析解在 0,1 区间的误差曲线，由图可知，在步长相同的情况下，梯形法的精确度略高于改进的欧拉法；改进的欧拉法和梯形法精确度都明显高于欧拉法例 2 用欧拉法、改进的欧拉法和 Runge-Kutta 法求解初值问题并比较三种方法的结果方程为 n 1的伯努利方程，可求得解析解为现用 MATLAB软件编程，用题目要求的方法求解，可得如下图示结果：图 2 （a）步长为 0.2 时 R-K 法和解析解比较图 2 （b） 步长为 0.2 时改进的 Euler 法和解析解比较图 2 （c）步长为 0.2 时欧拉法和解析解比较上图 2（

15、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 时，用欧拉法

16、、改进的欧拉法，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

17、,x,gy, figure（3）;plot（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:y（i+1）=y（i）+h*feval（f,x（i）,y（i）;end改进的欧拉法程序function gy=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 法程序function Ry=RK（f,

18、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 应用实例: 耐用消费新产品的销售规律模型一种新产品进入市场以后， 常常会经历销售量首先慢慢增加然后逐渐慢慢下降的一个过程，由此给出的时间销售坐标系下的曲线称为产品的生命曲线 , 它的形状是钟形的. 不过对于较耐用的消费品，情况

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

20、出了科学解释 .3.1.1 模型的建立将顾客获得信息大致可分成两类，其一称之为“搜集型” ，源自于产告说明、广告，“创新型”顾客获得这类消息后就可作出其是否购买；第二类称之为“体验型” ，即消费者使用之后会获得真实体验，常常以口头的形式散播， “模仿型”类顾客会在获得这种信息之后才能决定是否购买 .设 K 是潜在用户的总数， K1 与K2 分别为 “创新型”与“模仿型”的人数 . 再设N t 是时刻 t 已经购买的商品顾客的人数，而N t 与 N 2 t 分别为“创新型”与“模仿型”的顾客的人数 , 再设A t 是时刻 t 中已获得“搜集型”信息的人数，由于该部分的信息可直接由外部获得， 或者

21、可从已获得该类信息的人群中获得， 因此类似巴斯模型，从而建立如下方程：dA tdtK1 A1 t 1 2A1 t , A 0 0, 1, 2 0 ，获得“搜集型”信息的“创新型”消费者决定其是否购买的行为，满足如下方程：对于“模仿型”的顾客，可从已经购买该产品“创新型”或者“模仿型”的顾客中获取信息，于是有在这里，忽略消费者购买产品后需一段短暂的使用后才会散播体验信息的滞后作用综上，斯蒂芬斯莫赛模型是一常微分方程组的初值问题：而N t N t N t 为时刻 t 购买该商品的总人数 .3.1.2 模型的求解对于斯蒂芬斯莫赛模型中N t 的解析解则不能求出，于是可以用四阶Runge Kutta

22、公式求得，且在它的精度要求达到很高情形下求出 N2 t .用MATLAB软件求解，相应的程序及结果如下function RK=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;function f=

23、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）图 4N t N t 与时间关系图1 , 2图 5 N t 与0,25 时间段关系图由此例可以看出， 微分方程数值解在实际生活有着广泛的应用， 而数值解法中的Ru

24、nge-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-33