关闭

关闭

关闭

封号提示

内容

首页 第15章 常微分方程的解法.pdf

第15章 常微分方程的解法.pdf

第15章 常微分方程的解法.pdf

上传者: Gingerjin 2012-07-26 评分 0 0 0 0 0 0 暂无简介 简介 举报

简介:本文档为《第15章 常微分方程的解法pdf》,可适用于工程科技领域,主题内容包含第十五章常微分方程的解法建立微分方程只是解决问题的第一步通常需要求出方程的解来说明实际现象并加以检验。如果能得到解析形式的解固然是便于分析和应用的但符等。

第十五章常微分方程的解法建立微分方程只是解决问题的第一步通常需要求出方程的解来说明实际现象并加以检验。如果能得到解析形式的解固然是便于分析和应用的但是我们知道只有线性常系数微分方程并且自由项是某些特殊类型的函数时才可以肯定得到这样的解而绝大多数变系数方程、非线性方程都是所谓“解不出来”的即使看起来非常简单的方程如xydxdy=于是对于用微分方程解决实际问题来说数值解法就是一个十分重要的手段。常微分方程的离散化下面主要讨论一阶常微分方程的初值问题其一般形式是==)(),(yaybxayxfdxdy()在下面的讨论中我们总假定函数),(yxf连续且关于y满足李普希兹(Lipschitz)条件即存在常数L使得|||),(),(|yyLyxfyxf这样由常微分方程理论知初值问题()的解必定存在唯一。所谓数值解法就是求问题()的解)(xy在若干点bxxxxaN=<<<<=L处的近似值),,,(NnynL=的方法),,,(NnynL=称为问题()的数值解nnnxxh=称为由nx到nx的步长。今后如无特别说明我们总取步长为常量h。建立数值解法首先要将微分方程离散化一般采用以下几种方法:(i)用差商近似导数若用向前差商hxyxynn)()(代替)('nxy代入()中的微分方程则得),,())(,()()(L=nxyxfhxyxynnnn化简得))(,()()(nnnnxyxhfxyxy如果用)(nxy的近似值ny代入上式右端所得结果作为)(nxy的近似值记为ny则有),,(),(L==nyxhfyynnnn()这样问题()的近似解可通过求解下述问题===)(),,(),(ayynyxhfyynnnnL()得到按式()由初值y可逐次算出L,,yy。式()是个离散化的问题称为差分方程初值问题。需要说明的是用不同的差商近似导数将得到不同的计算公式。(ii)用数值积分方法将问题()的解表成积分形式用数值积分方法离散化。例如对微分方程两端积分得==),,())(,()()(nnxxnnndxxyxfxyxyL()右边的积分用矩形公式或梯形公式计算。(iii)Taylor多项式近似将函数)(xy在nx处展开取一次Taylor多项式近似则得))(,()()(')()(nnnnnnxyxhfxyxhyxyxy=再将)(nxy的近似值ny代入上式右端所得结果作为)(nxy的近似值ny得到离散化的计算公式),(nnnnyxhfyy=以上三种方法都是将微分方程离散化的常用方法每一类方法又可导出不同形式的计算公式。其中的Taylor展开法不仅可以得到求数值解的公式而且容易估计截断误差。欧拉(Euler)方法Euler方法Euler方法就是用差分方程初值问题()的解来近似微分方程初值问题()的解即由公式()依次算出)(nxy的近似值),,(L=nyn。这组公式求问题()的数值解称为向前Euler公式。如果在微分方程离散化时用向后差商代替导数即hxyxyxynnn)()()('则得计算公式===)(),,(),(ayynyxhfyynnnnL()用这组公式求问题()的数值解称为向后Euler公式。向后Euler法与Euler法形式上相似但实际计算时却复杂得多。向前Euler公式是显式的可直接求解。向后Euler公式的右端含有ny因此是隐式公式一般要用迭代法求解迭代公式通常为===),,,(),(),()()()(Lkyxhfyyyxhfyyknnnknnnnn()Euler方法的误差估计对于向前Euler公式()我们看到当L,,=n时公式右端的ny都是近似的所以用它计算的ny会有累积误差分析累积误差比较复杂这里先讨论比较简单的所谓局部截断误差。假定用()式时右端的ny没有误差即)(nnxyy=那么由此算出))(,()(nnnnxyxhfxyy=()局部截断误差指的是按()式计算由nx到nx这一步的计算值ny与精确值)(nxy之差)(nnyxy。为了估计它由Taylor展开得到的精确值)(nxy是)()('')(')()(hOxyhxhyxyxynnnn=()()、()两式相减(注意到),('yxfy=)得)()()('')(hOhOxyhyxynnn=()即局部截断误差是h阶的而数值算法的精度定义为:若一种算法的局部截断误差为)(phO则称该算法具有p阶精度。显然p越大方法的精度越高。式()说明向前Euler方法是一阶方法因此它的精度不高。改进的Euler方法梯形公式利用数值积分方法将微分方程离散化时若用梯形公式计算式()中之右端积分即))(,())(,())(,(nnnnxxxyxfxyxfhdxxyxfnn并用,nnyy代替)(),(nnxyxy则得计算公式),(),(=nnnnnnyxfyxfhyy这就是求解初值问题()的梯形公式。直观上容易看出用梯形公式计算数值积分要比矩形公式好。梯形公式为二阶方法。梯形公式也是隐式格式一般需用迭代法求解迭代公式为===),,,(),(),(),()()()(Lkyxfyxfhyyyxhfyyknnnnnknnnnn()由于函数),(yxf关于y满足Lipschitz条件容易看出||||)()()()(knknknknyyhLyy其中L为Lipschitz常数。因此当<<hL时迭代收敛。但这样做计算量较大。如果实际计算时精度要求不太高用公式()求解时每步可以只迭代一次由此导出一种新的方法改进Euler法。改进Euler法按式()计算问题()的数值解时如果每步只迭代一次相当于将Euler公式与梯形公式结合使用:先用Euler公式求ny的一个初步近似值ny称为预测值然后用梯形公式校正求得近似值ny即==校正预测),(),(),(nnnnnnnnnnyxfyxfhyyyxhfyy()式()称为由Euler公式和梯形公式得到的预测校正系统也叫改进Euler法。为便于编制程序上机式()常改写成===)(),(),(qpnpnnqnnnpyyyyhxhfyyyxhfyy()改进Euler法是二阶方法。龙格库塔(RungeKutta)方法回到Euler方法的基本思想用差商代替导数上来。实际上按照微分中值定理应有),(')()(<<=θθhxyhxyxynnn注意到方程),('yxfy=就有))(,()()(hxyhxhfxyxynnnnθθ=()不妨记))(,(hxyhxfKnnθθ=称为区间,nnxx上的平均斜率。可见给出一种斜率K()式就对应地导出一种算法。向前Euler公式简单地取),(nnyxf为K精度自然很低。改进的Euler公式可理解为K取),(nnyxf),(nnyxf的平均值其中),(nnnnyxhfyy=这种处理提高了精度。如上分析启示我们在区间,nnxx内多取几个点将它们的斜率加权平均作为K就有可能构造出精度更高的计算公式。这就是龙格库塔方法的基本思想。首先不妨在区间,nnxx内仍取个点仿照()式用以下形式试一下<<===,),,(),()(βαβαλλhkyhxfkyxfkkkhyynnnnnn()其中βαλλ,,,为待定系数看看如何确定它们使()式的精度尽量高。为此我们分析局部截断误差)(nnyxy因为)(nnxyy=所以()可以化为=====)())(,())(,())(,())(,()('))(,()()(hOxyxfhkxyxhfxyxfhkxyhxfkxyxyxfkkkhxyynnynnxnnnnnnnnnβαβαλλ()其中k在点))(,(nnxyx作了Taylor展开。()式又可表为)()()(')()(hOfffhxhyxyyyxnnn=αβαλλλ注意到)()('')(')()(hOxyhxhyxyxynnnn=中fy='yxfffy=''可见为使误差)()(hOyxynn=只须令=λλ=αλ=αβ()待定系数满足()的()式称为阶龙格库塔公式。由于()式有个未知数而只有个方程所以解不唯一。不难发现若令==λλ==βα即为改进的Euler公式。可以证明在,nnxx内只取点的龙格库塔公式精度最高为阶。阶龙格库塔公式要进一步提高精度必须取更多的点如取点构造如下形式的公式:=====),(),(),(),()(hkhkhkyhxfkhkhkyhxfkhkyhxfkyxfkkkkkhyynnnnnnnnnnβββαββαβαλλλλ()其中待定系数iiiβαλ,,共个经过与推导阶龙格库塔公式类似、但更复杂的计算得到使局部误差)()(hOyxynn=的个方程。取既满足这些方程、又较简单的一组iiiβαλ,,可得=====),(),(),(),()(hkyhxfkhkyhxfkhkyhxfkyxfkkkkkhynnnnnnnnn()这就是常用的阶龙格库塔方法(简称RK方法)。线性多步法以上所介绍的各种数值解法都是单步法这是因为它们在计算ny时都只用到前一步的值ny单步法的一般形式是),,,(),,(==NnhyxhyynnnnLϕ()其中),,(hyxϕ称为增量函数例如Euler方法的增量函数为),(yxf改进Euler法的增量函数为)),(,(),(),,(yxhfyhxfyxfhyx=ϕ如何通过较多地利用前面的已知信息如rnnnyyy,,,L来构造高精度的算法计算ny这就是多步法的基本思想。经常使用的是线性多步法。让我们试验一下=r即利用,nnyy计算ny的情况。从用数值积分方法离散化方程的()式=))(,()()(nnxxnndxxyxfxyxy出发记nnnfyxf=),(),(=nnnfyxf式中被积函数))(,(xyxf用二节点),(nnfx),(nnfx的插值公式得到(因)nxx所以是外插。)()())(,(==nnnnnnnnnnnnfxxfxxhxxxxfxxxxfxyxf()此式在区间,nnxx上积分可得))(,(=nnxxfhfhdxxyxfnn于是得到)(=nnnnffhyy()注意到插值公式()的误差项含因子))((nnxxxx在区间,nnxx上积分后得出h故公式()的局部截断误差为)(hO精度比向前Euler公式提高阶。若取L,,=r可以用类似的方法推导公式如对于=r有)(=nnnnnnffffhyy()其局部截断误差为)(hO。如果将上面代替被积函数))(,(xyxf用的插值公式由外插改为内插可进一步减小误差。内插法用的是,,,rnnnyyyL取=r时得到的是梯形公式取=r时可得)(=nnnnnnffffhyy()与()式相比虽然其局部截断误差仍为)(hO但因它的各项系数(绝对值)大为减小误差还是小了。当然()式右端的nf未知需要如同向后Euler公式一样用迭代或校正的办法处理。一阶微分方程组与高阶微分方程的数值解法一阶微分方程组的数值解法设有一阶微分方程组的初值问题==)(),,,,('iimiiyayyyyxfyL),,,(miL=()若记Tmyyyy),,,(L=Tmyyyy),,,(L=Tmffff),,,(L=则初值问题()可写成如下向量形式==)(),('yayyxfy()如果向量函数),(yxf在区域D:mRybxa,连续且关于y满足Lipschitz条件即存在>L使得对,baxmRyy,都有),(),(yyLyxfyxf那么问题()在,ba上存在唯一解)(xyy=。问题()与()形式上完全相同故对初值问题()所建立的各种数值解法可全部用于求解问题()。高阶微分方程的数值解法高阶微分方程的初值问题可以通过变量代换化为一阶微分方程组初值问题。设有m阶常微分方程初值问题====)()()()()()(,,)(',)(),,',,(mmmmyayyayyaybxayyyxfyLL()引入新变量)(,,',===mmyyyyyyL问题()就化为一阶微分方程初值问题========)()()()(),,,(')(')(')('mmmmmmmmyayyyxfyyayyyyayyyyayyyLMM()然后用中的数值方法求解问题()就可以得到问题()的数值解。最后需要指出的是在化学工程及自动控制等领域中所涉及的常微分方程组初值问题常常是所谓的“刚性”问题。具体地说对一阶线性微分方程组)(xAydxdyΦ=()其中ARym,,Φ为m阶方阵。若矩阵A的特征值),,,(miiL=λ满足关系),,,(RemiiL=<λ|Re|min|Re|maximiimiλλ>>则称方程组()为刚性方程组或Stiff方程组称数|Re|min|Re|maximiimisλλ=为刚性比。对刚性方程组用前面所介绍的方法求解都会遇到本质上的困难这是由数值方法本身的稳定性限制所决定的。理论上的分析表明求解刚性问题所选用的数值方法最好是对步长h不作任何限制。Matlab解法Matlab数值解非刚性常微分方程的解法Matlab的工具箱提供了几个解非刚性常微分方程的功能函数如odeodeode其中ode采用四五阶RK方法是解非刚性常微分方程的首选方法ode采用二三阶RK方法ode采用的是多步法效率一般比ode高。Matlab的工具箱中没有Euler方法的功能函数。(I)对简单的一阶方程的初值问题==)(),('yxyyxfy改进的Euler公式为===)(),(),(qpnpnnqnnnpyyyyhxhfyyyxhfyy我们自己编写改进的Euler方法函数eulerprom如下:functionx,y=eulerpro(fun,x,xfinal,y,n)ifnargin<,n=endh=(xfinalx)nx()=xy()=yfori=:nx(i)=x(i)hy=y(i)h*feval(fun,x(i),y(i))y=y(i)h*feval(fun,x(i),y)y(i)=(yy)end例用改进的Euler方法求解xxyy'=)(x)(=y解编写函数文件dotym如下:functionf=doty(x,y)f=*y*x^*x在Matlab命令窗口输入:x,y=eulerpro('doty',,,,)即可求得数值解。(II)odeodeode的使用Matlab的函数形式是t,y=solver('F',tspany)这里solver为odeodeode输入参数F是用M文件定义的微分方程),('yxfy=右端的函数。tspan=t,tfinal是求解区间y是初值。例用RK方法求解xxyy'=)(x)(=y解同样地编写函数文件dotym如下:functionf=doty(x,y)f=*y*x^*x在Matlab命令窗口输入:x,y=ode('doty',,,)即可求得数值解。刚性常微分方程的解法Matlab的工具箱提供了几个解刚性常微分方程的功能函数如odesodesodetodetb这些函数的使用同上述非刚性微分方程的功能函数。高阶微分方程),,',,()()(=nnyyytfyL解法例考虑初值问题)('')(')(''''''====yyyyyyy解(i)如果设'',',yyyyyy===那么======)(')(')('yyyyyyyyyyy初值问题可以写成)(),,('YYYtFY==的形式其中yyyY=。(ii)把一阶方程组写成接受两个参数t和y返回一个列向量的M文件Fm:functiondy=F(t,y)dy=y()y()*y()y()*y()注意:尽管不一定用到参数t和yM文件必须接受此两参数。这里向量dy必须是列向量。(iii)用Matlab解决此问题的函数形式为T,Y=solver('F',tspan,y)这里solver为ode、ode、ode输入参数F是用M文件定义的常微分方程组tspan=ttfinal是求解区间y是初值列向量。在Matlab命令窗口输入T,Y=ode('F',,)就得到上述常微分方程的数值解。这里Y和时刻T是一一对应的Y(:,)是初值问题的解Y(:,)是解的导数Y(:,)是解的二阶导数。例求vanderPol方程')(''=yyyyμ的数值解这里>μ是一参数。解(i)化成常微分方程组。设',yyyy==则有==)(''yyyyyyμ(ii)书写M文件(对于=μ)vdpm:functiondy=vdp(t,y)dy=y()(y()^)*y()y()(iii)调用Matlab函数。对于初值)(',)(==yy解为T,Y=ode('vdp',,)(iv)观察结果。利用图形输出解的结果:plot(T,Y(:,),'',T,Y(:,),'')title('SolutionofvanderPolEquation,mu=')xlabel('timet')ylabel('solutiony')legend('y','y')SolutionofvanderPolEquation,mu=timetsolutionyyy例vanderPol方程=μ(刚性)解(i)书写M文件vdpm:functiondy=vdp(t,y)dy=y()*(y()^)*y()y()(ii)观察结果t,y=odes('vdp',,)plot(t,y(:,),'o')title('SolutionofvanderPolEquation,mu=')xlabel('timet')ylabel('solutiony(:,)')常微分方程的解析解在Matlab中符号运算工具箱提供了功能强大的求解常微分方程的符号运算命令dsolve。常微分方程在Matlab中按如下规定重新表达:符号D表示对变量的求导。Dy表示对变量y求一阶导数当需要求变量的n阶导数时用Dn表示Dy表示对变量y求阶导数。由此常微分方程yyy='''在Matlab中将写成'Dy*Dy=y'。求解常微分方程的通解无初边值条件的常微分方程的解就是该方程的通解。其使用格式为:dsolve('diffequation')dsolve('diffequation''var')式中diffequation为待解的常微分方程第种格式将以变量t为自变量进行求解第种格式则需定义自变量var。例试解常微分方程')(=yyxyx解编写程序如下:symsxydiffequ='x^y(x*y)*Dy='dsolve(diffequ,'x')求解常微分方程的初边值问题求解带有初边值条件的常微分方程的使用格式为:dsolve('diffequation''conditioncondition…''var')其中conditioncondition…即为微分方程的初边值条件。例试求微分方程xyy=''''')('',)(',)(===yyy的解。解编写程序如下:y=dsolve('DyDy=x','y()=,Dy()=,Dy()=','x')求解常微分方程组求解常微分方程组的命令格式为:dsolve('diffequdiffequ…''var')dsolve('diffequdiffequ…''conditioncondition…''var')第种格式用于求解方程组的通解第种格式可以加上初边值条件用于具体求解。例试求常微分方程组:==xfgxgfcos''sin''的通解和在初边值条件为)(,)(,)('===gff的解。解编写程序如下:clc,clearequ='Df*g=sin(x)'equ='DgDf=cos(x)'generalf,generalg=dsolve(equ,equ,'x')f,g=dsolve(equ,equ,'Df()=,f()=,g()=','x')求解线性常微分方程组(i)一阶齐次线性微分方程组AXX='=nxxXM=nnnnaaaaALMMML这里的’表示对t求导数。Ate是它的基解矩阵。AXX=')(XtX=的解为)()(XetXttA=。例试解初值问题XX='=)(X解编写程序如下:symsta=,,,,,,x=x=expm(a*t)*x(ii)非齐次线性方程组由参数变易法可求得初值问题)('tfAXX=)(XtX=的解为=ttstAttAdssfeXetX)()()()(例试解初值问题=teXXtcos'=)(X。解编写程序如下:clc,clearsymstsa=,,,,,,ft=exp(t)*cos(*t)x=x=expm(a*t)*xint(expm(a*(ts))*subs(ft,s),s,,t)x=simple(x)习题十五用欧拉方法和龙格库塔方法求微分方程数值解画出解的图形对结果进行分析比较。(i)πππ',,)('''===yyynxxyyx(Bessel方程令)=n精确解xxyπsin=。(ii))(',)(,cos''===yyxyy幂级数解L=!!!!xxxxy一只小船渡过宽为d的河流目标是起点A正对着的另一岸B点。已知河水流速v与船在静水中的速度v之比为k。(i)建立小船航线的方程求其解析解。(ii)设=dm=vms=vms用数值解法求渡河所需时间、任意时刻小船的位置及航行曲线作图并与解析解比较。

用户评论(0)

0/200

精彩专题

上传我的资料

每篇奖励 +2积分

资料评价:

/13
0下载券 下载 加入VIP, 送下载券

意见
反馈

立即扫码关注

爱问共享资料微信公众号

返回
顶部