常微分方程模型的解

第五章微分方程模型,鲁东大学数学与信息学院魏建新wjx0426,第五章微分方程模型,5.1小序,引例5.2微分方程的解,5.1小序引例,在许多实际问题中,当直接导出变量之间的函数关系较为困难,但导出包含未知函数的导数或微分的关系式较为容易时,可用建立微分方程模型的方法来研究该问题.,常用建模方法,已知定律,模拟近似法,微元法,例1一个较热的物体置于室温为180c的房间内,该物体最初的温度是600c,3分钟以后降到500c.想知道它的温度降到300c需要多少间?10分钟以后它的温度是多少?,应用已知物理定律,牛顿冷却(加热)定律:将温度为T的物体放入处于常温m的介质中时,T的变化速率正比于T与周围介质的温度差.,分析:假设房间足够大,放入温度较低或较高的物体时,室内温度基本不受影响,即室温分布均衡,保持为m,采用牛顿冷却定律是一个相当好的近似.,建立模型:设物体在冷却过程中的温度为T(t),t0,常常事半功倍.,“T的变化速率正比于T与周围介质的温度差”,翻译为,数学语言,建立微分方程,其中参数k0，m=18.求得一般解为,jwfc.m,代入条件，求得,该物体温度降至300c需要13.8274分钟.,jwfun.m,jwkz.m,jwds.m,jwfun2.m,5.2.1微分方程的解析解5.2.2微分方程的数值解,5.2微分方程的解,在Matlab中,符号运算工具箱提供了功能强大的求解常微分方程的符号运算命令dsolve.格式：r=dsolve(eq1,eq2,.,cond1,cond2,.,v)r=dsolve(eq1,eq2,.,cond1,cond2,.,v)r=dsolve(方程1,方程n,初始条件,自变量)记号:在表达微分方程时，用字母D表示求微分，D2、D3等表示求高阶微分.任何D后所跟的字母为因变量，自变量可以指定或由系统规则选定为确省.,5.2.1微分方程的解析解,几个求解例子,结果：u=tg(t-c),ToMatlabff1,求解常微分方程的通解,解：输入命令:y=dsolve(D2y+4*Dy+29*y=0,y(0)=0,Dy(0)=15,x),结果为：y=3e-2xsin（5x）,ToMatlab(ff2）,求解常微分方程的初边值问题,解输入命令：x,y,z=dsolve(Dx=2*x-3*y+3*z,Dy=4*x-5*y+3*z,Dz=4*x-4*y+2*z,t)；x=simple(x)%将x化简y=simple(y)z=simple(z),ToMatlab(ff3),求解常微分方程组,例4,Tofeichangxishu1,Tofeichangxishu,Towujiexijie,例5,例6,y=dsolve(D2y+x*(y2-1)*Dy+y,x),练习,y=dsolve(D3y-D2y=x,y(1)=8,Dy(1)=7,D2y(2)=4,x),1求解,2.求解微分方程组的通解,及在初值条件下的解。,equ1=D2f+3*g=sin(x);equ2=Dg+Df=cos(x);general_f,general_g=dsolve(equ1,equ2,x)f,g=dsolve(equ1,equ2,Df(2)=0,f(3)=3,g(5)=1,x),5.2.2微分方程的数值解,建立微分方程只是解决问题的第一步,通常需要求出方程的解来说明实际现象,并加以检验.如果能得到解析形式的解固然是便于分析和应用的.但是我们知道,只有线性常系数微分方程,并且自由项是某些特殊类型的函数时,才可以肯定得到这样的解,而绝大多数变系数方程、非线性方程都是所谓“解不出来”的,即使看起来非常简单的方程,于是对于用微分方程解决实际问题来说,数值解法就是一个十分重要的手段.,（一）常微分方程数值解的定义,在生产和科研中所处理的微分方程往往很复杂且大多得不出一般解。而在实际上对初值问题，一般是要求得到解在若干个点上满足规定精确度的近似值，或者得到一个满足精确度要求的便于计算的表达式。,因此，研究常微分方程的数值解法是十分必要。,（二）建立数值解法的一些途径,1.1、用差商代替导数,若步长h较小，则有,故有公式：,此即欧拉法（显式Euler法）。,1.2、使用数值积分,对方程y=f(x,y),两边由xi到xi+1积分，并利用梯形公式，有：,实际应用时，与欧拉公式结合使用：,此即改进的欧拉法（梯形方法）。,故有公式：,1.3.预估-校正Euler法,综合使用前两种方法：用1去求初步的近似值，用之再代入2格式注：1.当一个数值公式的截断误差可表示为O(hk+1)时(k为正整数，h为步长),称它是一个k阶公式.2.欧拉法是一阶公式,改进的欧拉法是二阶公式.,2龙格库塔（RungeKutta）方法,2.1基本思想回到Euler方法的基本思想用差商代替导数上来.实际上,按照微分中值定理应有,注意到方程，有,不妨记,称之为的斜率.,如果在内多取几个点,将它们的斜率加权平均作为就有可能构造出精度更高的计算公式.这就是龙格库塔方法的基本思想。,2.2龙格-库塔方法一般形式,L称为“级”,其精度P称为“阶”,上式称为L级P阶龙格-库塔公式.,2.3经典4级4阶龙格-库塔方法,上面讨论的是一阶微分方程初值问题的数值解，这些方法可以直接推广到常微分方程组和高阶微分方程初值问题上.,（三）一阶微分方程组和高阶微分方程初值问题的数值解,设有一阶微分方程组的初值问题,若记,则上初值问题可写成如下向量形式,故对前面一阶微分方程初值问题所建立的各种数值解法可用于求解一阶微分方程组初值问题.,高阶微分方程的初值问题可以通过变量代换化为一阶微分方程组初值问题。,设有m阶常微分方程初值问题,引入新变量,化为一阶微分方程初值问题,(四)刚性常微分方程（组）,（五）用Matlab软件求常微分方程的数值解,t，x=solver（f,ts,x0,options）,常微分方程求解(不能确定是否是刚性常微分方程时，首先用命令ode45，然后用命令ode15S.),例题1求初值问题的数值解.,解：显然，解析解存在，为下面，求其数值解.,1建立函数文件fun1.mfunctiondy=fun1(x,y)dy=2*x;,2主程序shuzhijie.mx0=-1,1;y0=4X,Y=ode45(fun1,x0,4);plot(X,Y,+),shuzhijie.m,fun1.m,解:1、建立m-文件rigid.m如下：functiondy=rigid(t,y)dy=y(2)*y(3);-y(1)*y(3);-0.51*y(1)*y(2),2、取t0=0，tf=12，输入命令：T,Y=ode45(rigid,012,011);plot(T,Y(:,1),-,T,Y(:,2),*,T,Y(:,3),+),ToMatlab(ff4),例2解微分方程组.,例3：考虑初值问题,(2)把一阶方程组写成接受两个参数和返回一个列向量的M文件djl.m：functiondy=djl(t,y);dy=y(2);y(3);3*y(3)+y(2)*y(1);(或者dy=zeros(3,1);dy(1)=y(2);dy(2)=y(3);dy(3)=3*y(3)+y(2)*y(1);)(3)取t0=0，tf=1，输入命令：T,Y=ode45(rigid,01,01-1);plot(T,Y(:,1),-,T,Y(:,2),*,T,Y(:,3),+),ToMatlab(ff5),例4：求vanderPol方程的数值解，其中.,（2）书写M文件（对于）vdp1.m:,functiondy=vdp1(t,y);dy=y(2);(1-y(1)2)*y(2)-y(1);,（3）调用Matlab函数。对于初值,，解为T,Y=ode45(vdp1,020,2;0);plot(T,Y(:,1),-,T,Y(:,2),-),ToMatlab(ff7),1、建立m-文件vdp1000.m如下：functiond