微分方程数值解讲解

1、微分方程数值解及其应用绪论自然界中的许多事物的运动和变化规律都可以用微分方程来描述，因此对工程和 科学技术中的实际问题的研究中 , 常常需要求解微分方程 但往往只有少数较简单和 典型的微分方程可求出其解析解，在大多数情况下 , 只能用近似法求解，数值解法是 一类重要的近似方法 本文主要讨论一阶常微分方程的初值问题的数值解法， 探讨这 些算法在处理来自生活实际问题中的应用，并结合 MATLA软B 件，动手编程予以解决 微分方程的初值问题 11.1 预备知识在对生活实际问题的研究中，通常需要考虑一阶微分方程的初值问题dydxf (x, y)1）y(x0) y0这里 f x,y 是矩形区域 R ：

2、x x0 a, y y0 b 上的连续函数对初值问题（ 1）需要考虑以下问题：方程是否一定有解呢 ?若有解，有多少个解呢?下面给出相关的概念与定理定义 1 Lipschitz 条件 1 2：矩形区域 R ： x x0 a, y y0 b上的连续函数 f x, y 若满足：存在常数 L 0 ，使得不等式 f x, y1 f x,y2 L y1 y 2 对所有x,y1 , x,y2 R都成立，则称 f x,y 在R上关于 y 满足 Lipschitz 条件.定理 1 解的存在唯一性定理 1 3 ：设 f 在区域 Dx, y a x b,y R 上连 续，关于 y满足Lipschitz 条件，则对任

3、意的 x0 a,b ,y0 R，常微分方程初值问题（1）当 x a,b 时存在唯一的连续解 y x 该定理保证若一个函数 f x,y 关于 y满足Lipschitz 条件，它所对应的微分方程 的初值问题就有唯一解 . 在解的存在唯一性得到保证的前提下，自然要考虑方程的求解问题求解微分方程虽然有多种解析方法， 但根据工程和科学实践问题所得到的微 分方程往往很复杂，在很多情况下不能或很难给出解析解，有时即使能求出形式解， 也往往因形式过于复杂或计算量太大而不实用， 因此从实际问题中归结出来的微分方 程主要依靠数值解法定义 2 微分方程数值解：对初值问题（ 1）寻求数值解就是寻求解 y x 在一 系

4、列离散节点x0 x1 x2xn xn 1上的近似解 y0,y1,y2, ,yn,yn 1, ，相邻两个节点的间距 hn xn 1 xn称为步长.在一 般情况下假定 hi h i 0,1, 为常数，这时节点为 xn x0 nh,n 0,1,2, 要求微分方程数值解，首先要建立数值算法，即对初值问题（1）中的方程离散化，建立求解数值解法的递推公式一类是计算 yn 1时只用到前一点的值 yn ，称为单 步法；另一类是用到 yn 1前面k点的值 yn,yn1, ,yn k1称为k步法.对初值问题（ 1）式的单步法可用一般形式表示为yn 1 yn h （xn,yn,yn 1,h），其中多元函数 与 f

5、x,y 有关，当 含有 yn 1 时，方法是隐式的；若 中不含 yn 1， 则为显式方法，所以显式单步法可表示为yn 1 yn h （xn,yn,h）.（2）设 y x 是初值问题（ 1）的准确解，称Tn 1yxn 1yxnh（xn,yxn,h）为显式单步法（ 2）的局部截断误差 . 若存在最大正整数 p, 使显式单步法（ 2）式的局部截断误差满足 Tn 1 y x h y x h x, y,h O hp 1 ，则称（ 2）式有 p阶精度 .1.2 几种常用的数值解法及其分析、比较1.2.1 欧拉法与后退欧拉法1）欧拉法：欧拉曾简单地用差分代替微分，即利用公式y(xn 1) y(xn)hy(x

6、n) f (xn,yn)将初值问题（ 1）离散化，则问题（ 1）可化为yn 1 yn h f （xn ,yn）, xn x0 n h ，（3）此方法称为欧拉法 .欧拉方法的几何意义在数值计算公式中体现了出来 . 在 xy平面上，一阶微分方 程的解 y y x 称作它的积分曲线 . 积分曲线上一点 x, y 的切线斜率等于函数 f x, y ，按函数 f x,y 在 xy平面上建立一个方向场，那么，积分曲线上每一点的切 线方向均与方向场在该点的方向相一致 .基于上述几何解释， 从初始点 P0（ x0, y0 ）出发，先依方向场在该点的方向上推进到 x x1上一点 P1，再从 P1依方向场的方向推

7、进到 x x2上一点 P2 ，循环前进便作出一 条折线 P0 P1P2 ，因此欧拉方法又称为折线法 .若初值 y0已知，则由（ 3）式可逐步算 出y1 y0 h f（x0, y0）,y2 y1 h f （x1, y1）,为了分析计算公式的精确度，通常可用泰勒展开将 y xn 1 在 xn 处展开，则有nxn,xn 1 .y xn 1yxnh yxnyxnhh2yn ,在 yn y xn 的前提下， f xn,ynfxn, y xny xn . 可得欧拉法（ 3）的误差为h2h2y xn 1 yn 1 h2 y n h2 y xn .容易看出，欧拉法（ 3）式具有一阶精度 .2）向后欧拉方法：如

8、果对微分方程（ 1）从 xn 到 xn 1积分，得4）xn 1y xn 1 y xnx f t,y t dx ，xn如果（ 4）式右端积分用右矩形公式 h f xn 1,y xn 1 近似，则得到另一个公式yn 1 yn hf xn1,yn 1 ，5）称为后退欧拉法 .值得一提的是 :后退欧拉法与欧拉公式有着本质的区别，后者是关于yn 1的直接计算公式，它是显式的，而( 5)式的右端含有关于 yn 1的表达式，它是隐式的 . 在利 用后退欧拉法时，我们通常利用迭代法求解，实质就是逐步显示化 . 具体迭代过程如 下：首先利用欧拉公式 yn(0)1 yn h f ( xn , yn )给出迭代初值

9、 yn(0)1，把它代入( 5)式的右端，使之转化为显式，直接计算得 y(1) yn h f (xn 1, y(0) ) . 如此反复进行，得y(nk11) yn h f (xn 1,y(nk)1) k 0,1, ，则得到后退欧拉法的迭代公式yn(0)1 yn h f (xn,yn)，yn(k11) yn h f (xn 1,yn(k)1)可以看出，后退欧拉法具有一阶精度，且计算比较麻烦 .1.2.2 梯形方法为得到比欧拉法精确度高的计算公式，在等式( 4)式右端积分中若用梯形求积公式近似，并用 yn 代替 y xn ， yn 1代替 y xn 1 ，则得6)yn 1 yn h2 f xn,y

10、n f xn 1,yn 1 ， 称其为梯形方法 .梯形方法与后退欧拉法一样， 都是隐式单步法， 可用迭代法求解， 其迭代公式为7)yn(0)1 yn h f (xn,yn)yn(k11) yn h2 f xn, yn f xn 1, yn(k)1为了分析梯形公式的收敛性，将( 6)与( 7)式相减，得yn 1 yn(k11) h2 f xn 1, yn 1 f xn 1, yn(k 1) ， k 0,1,2,因为 f x,y 满足 Lipschitz 条件，于是有yn 1 yn 1h2L yn 1 yn 1 ，其中 L 为 f x,y关 于 y 的 Lipschitz 常 数 . 如 果 选

11、取 h 充 分 小 ， 使 得 hL 1 ， 则 当 k时 有2y(nk11) yn 1，这说明迭代过程( 7)式是收敛的 4 .容易推导得出梯形法( 7)式是二阶方法.经分析，梯形方法虽然提高了精度，但是以增加计算量为代价的 . 从上述的迭代公式可以看出，每迭代一次都要重新计算 f x, y 的值，而且迭代又要进行若干次， 计算相当的复杂 . 为此，有没有比较简便的计算方法呢？下面给出改进的欧拉方法.1.2.3 改进的欧拉方法由前面的讨论可知， 梯形法计算相对复杂， 现对上面的梯形法进行简化， 具体方 法是只计算一两次就转入下一步的计算，先用欧拉公式（ 3）求得一个初步的近似解 yn 1，称

12、为预测值，再利用公式（ 6）把它校正一次，这样建立的预测 -校正系统通常 称为改进的欧拉公式 . 具体公式如下hyn 1 ynf xn,yn f xn 1,yn hf xn,yn（8）2改进的欧拉法与梯形法一样，是二阶方法 .1.2.4 Runge-Kutta 方法由前面讨论可知，从（ 4）式xn 1y xn 1 y xn x f x, y x dxxn可以看出， 若要使公式阶数提高， 就必须使右端积分的数值求积公式精度提高， 它必 然要增加求积积点，为此将（ 4）式的右端用求积公式表示为xn 1rx f x,y x dx h ci f xn ih,y xn ih ，（9）xni 1 i n

13、i n i一般来说，点数 r 越多，精度越高，上式右端相当于增量函数 x,y,h ，为得到便于 计算的显式方法，将公式（ 9）表示为：yn 1 yn h xn,yn,h ,（10）其中rxn,yn,hciKii1K1 f xn,yn（11）i1Ki f xn ih,yn h ijK j ,i 2, rj1这里ci, i, ij 均为常数. ci为加权因子， Ki为第i段斜率，共有 r 段.我们把（ 10）和（11）称为 r 级显式 Runge-Kutta 法，简称为 R-K 方法 . 下面给出其中最经典最常用 的一个公式：hyn 1 ynK1 2K2 2K3 K 4 ,6K1xn,ynhhK2

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

15、y 的值，计算量较大（对于复杂的 f x,y 而言）2 数值方法的应用实例 5-9y 10y例1 对于初值问题 ，分别用欧拉法、 改进的欧拉法，梯形法求 y 1 的 y 0 1近似值.解：易得该方程的解析解 y x e 10x ， y 1 4.5400e-005，为比较，将按不同 数值计算方法所得结果列表如下：h欧拉法改进的欧拉法梯形法0.2-1110.109.7656E-0045.4994E-0050.012.6561 E-0054.6223 E-0054.5026 E -0050.0014.3717 E-0054.5408 E-0054.5396 E -0050.00014.5173 E-

16、0054.5400 E-0054.5400 E -005表 1 三种不同方法的数值结果0.020.0150.010.00500.2 0.4欧拉法误差 改进欧拉法误差 梯形法误差0.6 0.8图 1 三种不同方法数值解与精确解的误差曲线从表 1 中可以看出：当 h 0.2时，三种方法均不稳定， 计算结果严重偏离精确值； h 0.1时，改进后的欧拉和梯形法均稳定，但欧拉法效果很差；当 h 0.01 时，三种 方法均稳定，但精确度有区别可以看出， h 越小，计算结果越好，要想计算结果充 分接近于解析解还须取较小的 h 值图 1 反映的步长 h 0.01时，三种数值方法的所得数值解与解析解在 0,1

17、区间的误差 曲线，由图可知，在步长相同的情况下，梯形法的精确度略高于改进的欧拉法；改进 的欧拉法和梯形法精确度都明显高于欧拉法例 2 用欧拉法、改进的欧拉法和 Runge-Kutta 法求解初值问题dy y 2x,x 0,1dx yy 0 1并比较三种方法的结果解：方程为 n 1 的 伯努利方程，可求得解析解为y e2x e 4x 2xe 4x现用 MATLAB 软件编程，用题目要求的方法求解，可得如下图示结果：1.81.61.41.20.20.40.60.8解析解 Runge-kutta 法0图 2 (a)步长为 0.2 时 R-K 法和解析解比较1.81.61.41.2解析解改 进 的 E

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

19、uler 法1.61.41.210 0.2 0.4 0.6 0.8 1图 3 (a) 步长为 0.1 时 Euler 法和解析解比较1.81.61.41.21解析解改 进 的 Euler 法0.20.4 0.6 0.8图 3 (b) 步长为 0.1 时改进的 Euler 法和解析解比较1.81.61.41.21解析解Runge-kutta 法00.20.4 0.6 0.8图 3 (c) 步长为 0.1 时 Runge-Kutta 法和解析解比较上图 3 (a)，(b)，(c)描述的是步长为 0.1 时，用欧拉法、改进的欧拉法， Runge-Kutta 法求解方程所得的数值解与解析解之间的对比图

20、由图可知，改进的欧拉法和 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, *)

21、欧拉法程序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:ny(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:nyp=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;endgy=y;Runge-Kutta法

22、程序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:nk1=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);10 k4=feval(f,x(i+1),y(i)+h*k3);y(i+1)=y(i)+h*(k1+2*k2+2*k3+k4)/6;endRy=y;3 微分方程数值解法在实际生活中的应用3.1 应用实例 : 耐用消费新产品的销售规律模型一种新产品进入市场以后， 常常会经

23、历销售量首先慢慢增加然后逐渐慢慢下降的 一个过程，由此给出的时间销售坐标系下的曲线称为产品的生命曲线 , 它的形状是 钟形的. 不过对于较耐用的消费品，情况会有所不一样，它的生命曲线会在开始有个 小小的高峰，之后是段比较平坦的曲线，先下降，之后再上升，然后达到高峰，因此 它是双峰型的曲线 .如何理解这种和传统产品的生命曲线理论相冲突的现象？澳大利亚学者斯蒂芬 斯与莫赛观察到的购买耐用消费产品的人大概可分为两类： 其一是非常善于接受新的 事物，称其为“创新型”消费者，他们会经常从产品广告，制造商给出产品的说明书 与商店样品中了解了产品功能与性能之后， 再决定其否购买； 其二是消费者相对保守， 称

24、其为“模仿型”消费者，他们往往会根据大部分已经购买产品的消费者实际使用的 经验而提供的信息来决定是否购买产品， 下面的例子经过分析， 建立相应的数学模型， 对这种现象给出了科学解释 .3.1.1 模型的建立将顾客获得信息大致可分成两类，其一称之为“搜集型” ，源自于产告说明、广 告，“创新型”顾客获得这类消息后就可作出其是否购买；第二类称之为“体验型” ， 即消费者使用之后会获得真实体验，常常以口头的形式散播， “模仿型”类顾客会在 获得这种信息之后才能决定是否购买 .设 K 是潜在用户的总数， K1与 K2分别为 “创新型”与“模仿型”的人数 .再 设 N t 是时刻 t 已经购买的商品顾客

25、的人数，而 N 1 t 与 N 2 t 分别为“创新型”与 “模仿型”的顾客的人数 ,再设 A1 t 是时刻 t中已获得“搜集型”信息的人数，由于 该部分的信息可直接由外部获得， 或者可从已获得该类信息的人群中获得， 因此类似 巴斯模型，从而建立如下方程：11dA1 tK1 A1 t1 2A1 t ,A 0 0, 1, 2 0 ，dt获得“搜集型”信息的“创新型”消费者决定其是否购买的行为，满足如下方程：dN1 tK1 N1 tN1 t ,N1 0 0, , 0dt对于“模仿型”的顾客，可从已经购买该产品“创新型”或者“模仿型”的顾客中获 取信息，于是有dN2 tdtK2N2 t N1 t N

26、2 t , 0在这里，忽略消费者购买产品后需一段短暂的使用后才会散播体验信息的滞后作用综上，斯蒂芬斯莫赛模型是一常微分方程组的初值问题：dN1 tdtK1 N1 tN1 t ,dN2 t2 K2 N2 t N1 t N2 t , dtN1 0 0, N2(0) 0而 N t N1 t N2 t 为时刻 t 购买该商品的总人数 .3.1.2 模型的求解对于斯蒂芬斯莫赛模型中 N2 t 的解析解则不能求出，于是可以用四阶Runge Kutta 公式求得，且在它的精度要求达到很高情形下求出 N2 t .用 MATLAB软件求解，相应的程序及结果如下function RK=RKFC(fc,a,b,h,

27、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;endRK=y;12function 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.