微分方程数值解.doc

-微分方程数值解及其应用绪论自然界中的许多事物的运动和变化规律都可以用微分方程来描述，因此对工程和科学技术中的实际问题的研究中, 常常需要求解微分方程但往往只有少数较简单和典型的微分方程可求出其解析解，在大多数情况下,只能用近似法求解，数值解法是一类重要的近似方法本文主要讨论一阶常微分方程的初值问题的数值解法，探讨这些算法在处理来自生活实际问题中的应用，并结合MATLAB软件，动手编程予以解决微分方程的初值问题1.1 预备知识在对生活实际问题的研究中，通常需要考虑一阶微分方程的初值问题 1这里是矩形区域：上的连续函数对初值问题1需要考虑以下问题：方程是否一定有解呢"假设有解，有多少个解呢"下面给出相关的概念与定理定义1条件：矩形区域：上的连续函数假设满足：存在常数，使得不等式对所有都成立，那么称在上关于满足条件.定理 1解的存在唯一性定理：设在区域上连续，关于满足条件，那么对任意的，常微分方程初值问题1当时存在唯一的连续解该定理保证假设一个函数关于满足条件，它所对应的微分方程的初值问题就有唯一解.在解的存在唯一性得到保证的前提下，自然要考虑方程的求. z.-解问题求解微分方程虽然有多种解析方法，但根据工程和科学实践问题所得到的微分方程往往很复杂，在很多情况下不能或很难给出解析解，有时即使能求出形式解，也往往因形式过于复杂或计算量太大而不实用，因此从实际问题中归结出来的微分方程主要依靠数值解法定义 2微分方程数值解：对初值问题1寻求数值解就是寻求解在一系列离散节点上的近似解，相邻两个节点的间距称为步长.在一般情况下假定为常数，这时节点为要求微分方程数值解，首先要建立数值算法，即对初值问题1中的方程离散化，建立求解数值解法的递推公式一类是计算时只用到前一点的值，称为单步法；另一类是用到前面点的值称为步法.对初值问题1式的单步法可用一般形式表示为，其中多元函数与有关，当含有时，方法是隐式的；假设中不含，那么为显式方法，所以显式单步法可表示为. 2设是初值问题1的准确解，称为显式单步法2的局部截断误差. 假设存在最大正整数,使显式单步法2式的局部截断误差满足，那么称2式有阶精度.1.2几种常用的数值解法及其分析、比拟欧拉法与后退欧拉法1欧拉法：欧拉曾简单地用差分代替微分，即利用公式将初值问题1离散化，那么问题1可化为， 3此方法称为欧拉法.欧拉方法的几何意义在数值计算公式中表达了出来.在平面上，一阶微分方程的解称作它的积分曲线.积分曲线上一点的切线斜率等于函数，按函数在平面上建立一个方向场，那么，积分曲线上每一点的切线方向均与方向场在该点的方向相一致.基于上述几何解释，从初始点出发，先依方向场在该点的方向上推进到上一点，再从依方向场的方向推进到上一点，循环前进便作出一条折线，因此欧拉方法又称为折线法.假设初值，那么由3式可逐步算出 为了分析计算公式的准确度，通常可用泰勒展开将在处展开，那么有在的前提下，可得欧拉法3的误差为容易看出，欧拉法3式具有一阶精度.2向后欧拉方法：如果对微分方程1从到积分，得， 4如果4式右端积分用右矩形公式近似，那么得到另一个公式 ， 5称为后退欧拉法.值得一提的是:后退欧拉法与欧拉公式有着本质的区别，后者是关于的直接计算公式，它是显式的，而5式的右端含有关于的表达式，它是隐式的.在利用后退欧拉法时，我们通常利用迭代法求解，实质就是逐步显示化.具体迭代过程如下：首先利用欧拉公式给出迭代初值，把它代入5式的右端，使之转化为显式，直接计算得.如此反复进展，得 ，那么得到后退欧拉法的迭代公式，可以看出，后退欧拉法具有一阶精度，且计算比拟麻烦.梯形方法为得到比欧拉法准确度高的计算公式，在等式4式右端积分中假设用梯形求积公式近似，并用代替，代替，那么得， 6称其为梯形方法.梯形方法与后退欧拉法一样，都是隐式单步法，可用迭代法求解，其迭代公式为. 7为了分析梯形公式的收敛性，将6与7式相减，得，因为满足条件，于是有，其中为关于的常数.如果选取充分小，使得，那么当时有，这说明迭代过程7式是收敛的4.容易推导得出梯形法7式是二阶方法.经分析，梯形方法虽然提高了精度，但是以增加计算量为代价的.从上述的迭代公式可以看出，每迭代一次都要重新计算的值，而且迭代又要进展假设干次，计算相当的复杂.为此，有没有比拟简便的计算方法呢？下面给出改良的欧拉方法.改良的欧拉方法由前面的讨论可知，梯形法计算相对复杂，现对上面的梯形法进展简化，具体方法是只计算一两次就转入下一步的计算，先用欧拉公式3求得一个初步的近似解，称为预测值，再利用公式6把它校正一次，这样建立的预测-校正系统通常称为改良的欧拉公式.具体公式如下 8改良的欧拉法与梯形法一样，是二阶方法.Runge-Kutta方法由前面讨论可知，从4式可以看出，假设要使公式阶数提高，就必须使右端积分的数值求积公式精度提高，它必然要增加求积积点，为此将4式的右端用求积公式表示为， 9一般来说，点数越多，精度越高，上式右端相当于增量函数，为得到便于计算的显式方法，将公式9表示为： 10其中 11 这里均为常数.为加权因子，为第段斜率，共有段.我们把10和11称为级显式Runge-Kutta法，简称为R-K方法.下面给出其中最经典最常用的一个公式： 12Runge-Kutta方法作为一种重要的单步方法，具有很高的实用价值，它关于初值是稳定的，其解连续地依赖于初值，是一类便于应用的单步法，为了计算，只用到前面一步的值即可，因此每步的步长可以独立取定.常用的Runge-Kutta方法精度较高，为了到达预定的精度，与欧拉方法与梯形法相比，步长可取得大些，求解区间上的总步数可以少些.但Runge-Kutta方法也有些缺点，比方四阶Runge-Kutta方法每算一步需要四次计算的值，计算量较大对于复杂的而言2数值方法的应用实例5-9例1对于初值问题，分别用欧拉法、改良的欧拉法，梯形法求的近似值.解：易得该方程的解析解，为比拟，将按不同数值计算方法所得结果列表如下：表 1三种不同方法的数值结果欧拉法改良的欧拉法梯形法0.2-1110.109.7656E-0045.4994E-0050.012.6561-0054.6223-0054.5026-0050.0014.3717-0054.5408-0054.5396-0050.00014.5173-0054.5400-0054.5400-005图 1三种不同方法数值解与准确解的误差曲线从表1中可以看出：当时，三种方法均不稳定，计算结果严重偏离准确值；时，改良后的欧拉和梯形法均稳定，但欧拉法效果很差；当时，三种方法均稳定，但准确度有区别可以看出，越小，计算结果越好，要想计算结果充分接近于解析解还须取较小的值图1反映的步长时，三种数值方法的所得数值解与解析解在区间的误差曲线，由图可知，在步长一样的情况下，梯形法的准确度略高于改良的欧拉法；改良的欧拉法和梯形法准确度都明显高于欧拉法例2用欧拉法、改良的欧拉法和Runge-Kutta法求解初值问题并比拟三种方法的结果解：方程为的伯努利方程，可求得解析解为现用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);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: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法程序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);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应用实例:耐用消费新产品的销售规律模型 一种新产品进入市场以后，常常会经历销售量首先慢慢增加然后逐渐慢慢下降的一个过程，由此给出的时间销售坐标系下的曲线称为产品的生命曲线,它的形状是钟形的.不过对于较耐用的消费品，情况会有所不一样，它的生命曲线会在开场有个小小的顶峰，之后是段比拟平坦的曲线，先下降，之后再上升，然后到达顶峰，因此它是双峰型的曲线. 如何理解这种和传统产品的生命曲线理论相冲突的现象？澳大利亚学者斯蒂芬斯与莫赛观察到的购置耐用消费产品的人大概可分为两类：其一是非常善于承受新的事物，称其为"创新型消费者，他们会经常从产品广告，制造商给出产品的说明书与商店样品中了解了产品功能与性能之后，再决定其否购置；其二是消费者相对保守，称其为"模仿型消费者，他们往往会根据大局部已经购置产品的消费者实际使用的经历而提供的信息来决定是否购置产品，下面的例子经过分析，建立相应的数学模型，对这种现象给出了科学解释.3.1.1 模型的建立 将顾客获得信息大致可分成两类，其一称之为"搜集型，源自于产告说明、广告，"创新型顾客获得这类消息后就可作出其是否购置；第二类称之为"体验型，即消费者使用之后会获得真实体验，常常以口头的形式散播，"模仿型类顾客会在获得这种信息之后才能决定是否购置. 设是潜在用户的总数，与分别为 "创新型与"模仿型的人数.再设是时刻已经购置的商品顾客的人数，而与分别为"创新型与"模仿型的顾客的人数,再设是时刻中已获得"搜集型信息的人数，由于该局部的信息可直接由外部获得，或者可从已获得该类信息的人群中获得，因此类似巴斯模型，从而建立如下方程：，获得"搜集型信息的"创新型消费者决定其是否购置的行为，满足如下方程：对于"模仿型的顾客，可从已经购置该产品"创新型或者"模仿型的顾客中获取信息，于是有在这里，忽略消费者购置产品后需一段短暂的使用后才会散播体验信息的滞后作用综上，斯蒂芬斯莫赛模型是一常微分方程组的初值问题：而为时刻购置该商品的总人数.3.1.2 模型的求解对于斯蒂芬斯莫赛模型中的解析解那么不能求出，于是可以用四阶公式求得，且在它的精度要求到达很高情形下求出.用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: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;endRK=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 与时间关系图图 5 与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 puting.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. z.