Matlab求解微分方程(组)及偏微分方程(组)教学教材

第四讲【2 】Matlab求解微分方程（组）理论介绍：Matlab求解微分方程（组）敕令求解实例：Matlab求解微分方程（组）实例现实运用问题经由过程数学建模所归纳得到的方程,绝大多半都是微分方程,真正能得到代数方程的机遇很少.另一方面,可以或许求解的微分方程也是十分有限的,特别是高阶方程和偏微分方程（组）.这就请求我们必须研讨微分方程（组）的解法：解析解法和数值解法.一．相干函数.敕令及简介1.在Matlab中,用大写字母D表示导数,Dy表示y关于自变量的一阶导数,D2y表示y关于自变量的二阶导数,依此类推.函数dsolve用来解决常微分方程（组）的求解问题,挪用格局为：X=dsolve(‘eqn1’,’eqn2’,…)函数dsolve用来解符号常微分方程.方程组,假如没有初始前提,则求出通解,假如有初始前提,则求出特解.留意,体系缺省的自变量为t2.函数dsolve求解的是常微分方程的准确解法,也称为常微分方程的符号解.但是,有大量的常微分方程固然从理论上讲,其解是消失的,但我们却无法求出其解析解,此时,我们须要追求方程的数值解,在求常微分方程数值解方面,MATLAB 具有丰硕的函数,我们将其统称为solver,其一般格局为：[T,Y]=solver(odefun,tspan,y0)解释：(1)solver为敕令ode45.ode23.ode113.ode15s.ode23s.ode23t.ode23tb.ode15i之一.(2)odefun 是显示微分方程'(,)y f t y =在积分区间tspan 0[,]f t t =上从0t 到f t 用初始前提y 求解.(3)假如要获得微分方程问题在其他指准时光点012,,,,ft t t t 上的解,则令tspan012[,,,]f t t t t =(请求是单调的).(4)因为没有一种算法可以有用的解决所有的ODE 问题,为此,Matlab 供给了多种求解器solver,对于不同的ODE 问题,采用不同的solver.表1 Matlab 中文本文件读写函数解释：ode23.ode45是极其常用的用来求解非刚性的标准情势的一阶微分方程（组）的初值问题的解的Matlab 常用程序,个中：ode23采用龙格-库塔2阶算法,用3阶公式作误差估量来调节步长,具有低等的精度.ode45则采用龙格-库塔4阶算法,用5阶公式作误差估量来调节步长,具有中等的精度.3．在matlab敕令窗口.程序或函数中创建局部函数时,可用内联函数inline,inli ne函数情势相当于编写M函数文件,但不需编写M-文件就可以描写出某种数学关系.挪用inline函数,只能由一个matlab表达式构成,并且只能返回一个变量,不许可[u,v]这种向量情势.因而,任何请求逻辑运算或乘法运算以求得最终成果的场合,都不能运用inline函数,inline函数的一般情势为：FunctionName=inline(‘函数内容’, ‘所有自变量列表’)例如：（求解F(x)=x^2*cos(a*x)-b ,a,b是标量;x是向量）在敕令窗口输入: Fofx=inline(‘x .^2*cos(a*x)-b’ , ‘x’,’a’,’b’);g= Fofx([pi/3 pi/3.5],4,1)体系输出为：g=-1.5483 -1.7259留意：因为运用内联对象函数inline不须要别的树立m文件,所有运用比较便利,别的在运用ode45函数的时刻,界说函数往往须要编辑一个m文件来单独界说,如许不便于治理文件,这里可以运用inline来界说函数.二．实例介绍1.几个可以直接用Matlab求微分方程准确解的实例例1 求解微分方程2 '2x y xy xe-+=程序：syms x y; y=dsolve(‘Dy+2*x*y=x*exp(-x^2)’,’x’)例 2求微分方程'0xxy y e+-=在初始前提(1)2y e=下的特解并画出解函数的图形.程序：syms x y; y=dsolve(‘x*Dy+y-exp(1)=0’,’y(1)=2*exp(1)’,’x’);ezplot(y)例 3求解微分方程组530tdx x y e dt dy x y dt ⎧++=⎪⎪⎨⎪--=⎪⎩在初始前提00|1,|0t t x y ====下的特解并画出解函数的图形.程序：syms x y t[x,y]=dsolve('Dx+5*x+y=exp(t)','Dy-x-3*y=0','x(0)=1','y(0)=0','t') simple(x); simple(y)ezplot(x,y,[0,1.3]);axis auto2.用ode23.ode45等求解非刚性标准情势的一阶微分方程（组）的初值问题的数值解（近似解）例4求解微分方程初值问题2222(0)1dy y x x dx y ⎧=-++⎪⎨⎪=⎩的数值解,求解规模为区间[0,0.5].程序：fun=inline('-2*y+2*x^2+2*x','x','y'); [x,y]=ode23(fun,[0,0.5],1); plot(x,y,'o-') 例5求解微分方程22'2(1)0,(0)1,(0)0d y dy y y y y dt dt μ--+===的解,并画出解的图形.剖析：这是一个二阶非线性方程,我们可以经由过程变换,将二阶方程化为一阶方程组求解.令12,,7dyx y x dt μ===,则121221212,(0)17(1),(0)0dx x x dtdx x x x x dt ⎧==⎪⎪⎨⎪=--=⎪⎩编写M-文件vdp.m function fy=vdp(t,x)fy=[x(2);7*(1-x(1)^2)*x(2)-x(1)]; end在Matlab 敕令窗口编写程序 y0=[1;0][t,x]=ode45(@vdp,[0,40],y0);或[t,x]=ode45('vdp',[0,40],y0); y=x(:,1);dy=x(:,2); plot(t,y,t,dy)演习与思虑：M-文件vdp.m 改写成inline 函数程序？ 3.用Euler 折线法求解Euler 折线法求解的根本思惟是将微分方程初值问题00(,)()dyf x y dx y x y ⎧=⎪⎨⎪=⎩化成一个代数(差分)方程,重要步骤是用差商()()y x h y x h +-替代微商dydx ,于是00()()(,())()k k k k y x h y x f x y x h y y x +-⎧=⎪⎨⎪=⎩记1,(),k k k k x x h y y x +=+=从而1(),k k y y x h +=+于是0011(),,0,1,2,,1(,).k k k k k k y y x x x h k n y y hf x y ++=⎧⎪=+=-⎨⎪=+⎩例6用Euler 折线法求解微分方程初值问题22(0)1dyx y dxy y ⎧=+⎪⎨⎪=⎩的数值解（步长h 取0.4）,求解规模为区间[0,2].剖析：本问题的差分方程为00110,1,0.4,0,1,2,,1(,).k k k k k k x y h x x h k n y y hf x y ++===⎧⎪=+=-⎨⎪=+⎩程序：>> clear >> f=sym('y+2*x/y^2'); >> a=0; >> b=2; >> h=0.4; >> n=(b-a)/h+1; >> x=0; >> y=1;>> szj=[x,y];%数值解 >> for i=1:n-1y=y+h*subs(f,{'x','y'},{x,y});%subs,调换函数 x=x+h; szj=[szj;x,y];end >>szj>> plot(szj(:,1),szj(:,2)) 解释：调换函数subs 例如：输入subs(a+b,a,4)意思就是把a 用4调换掉落,返回4+b,也可以调换多个变量,例如：subs(cos(a)+sin(b),{a,b},[sym('alpha'),2])分离用字符alpha 调换a 和2调换b,返回 cos(alpha)+sin(2)特别解释：本问题可进一步运用四阶Runge-Kutta 法求解,Euler 折线法现实上就是一阶Runge-Kutta 法,Runge-Kutta 法的迭代公式为001112341213243(),,(22),6(,),0,1,2,,1(,),22(,),22(,).k k k k k k k k k k k k y y x x x h h y y L L L L L f x y k n h h L f x y L h h L f x y L L f x h y hL ++=⎧⎪=+⎪⎪=++++⎪⎪=⎪=-⎨⎪=++⎪⎪⎪=++⎪⎪=++⎩响应的Matlab 程序为：>> clear >> f=sym('y+2*x/y^2'); >> a=0; >> b=2; >> h=0.4; >> n=(b-a)/h+1;>> x=0; >> y=1;>> szj=[x,y];%数值解 >> for i=1:n-1l1=subs(f,{'x','y'},{x,y});调换函数 l2=subs(f,{'x','y'},{x+h/2,y+l1*h/2}); l3=subs(f,{'x','y'},{x+h/2,y+l2*h/2}); l4=subs(f,{'x','y'},{x+h,y+l3*h}); y=y+h*(l1+2*l2+2*l3+l4)/6; x=x+h; szj=[szj;x,y]; end >>szj>> plot(szj(:,1),szj(:,2))演习与思虑：(1)ode45求解问题并比较差异.(2)运用Matlab 求微分方程(4)(3)''20y y y -+=的解. (3)求解微分方程''2',2(1)0,030,(0)1,(0)0y y y y x y y --+=≤≤==的特解.(4)运用Matlab 求微分方程初值问题2''''00(1)2,|1,|3x x x y xy y y ==+===的解.提示：尽可能多的斟酌解法三．微分方程转换为一阶显式微分方程组Matlab 微分方程解算器只能求解标准情势的一阶显式微分方程（组）问题,是以在运用ODE 解算器之前,我们须要做的第一步,也是最重要的一步就是借助状况变量将微分方程（组）化成Matlab 可接收的标准情势.当然,假如ODEs 由一个或多个高阶微分方程给出,则我们应先将它变换成一阶显式常微分方程组.下面我们以两个高阶微分方程组构成的ODEs 为例介绍若何将它变换成一个一阶显式微分方程组.Step1将微分方程的最高阶变量移到等式左边,其它移到右边,并按阶次从低到高分列.情势为：()'''(1)'''(1)()'''(1)'''(1)(,,,,,,,,,,)(,,,,,,,,,,)m m n n m n x f t x x x x y y y y y g t x x x xy y y y ----⎧=⎨=⎩ Step 2 为每一阶微分式选择状况变量,最高阶除外'''(1)123'''(1)123,,,,,,,,,m m n m m m m n x x x x x x x x x y x y x y x y --++++========留意：ODEs 中所有是因变量的最高阶次之和就是须要的状况变量的个数,最高阶的微分式不须要给它状况变量.Step 3 依据选用的状况变量,写出所有状况变量的一阶微分表达式''''122334123''12123,,,,(,,,,,),,(,,,,,)m m n m m m n m n x x x x x x x f t x x x x x x x g t x x x x +++++======演习与思虑：(1)求解微分方程组**'''3312*'''3312()()22x x x y x r r y y y x y r r μμμμμμ⎧+-=+--⎪⎪⎨⎪=+--⎪⎩个中2r =1r =*1,μμ=-1/82.45,μ=(0) 1.2,x =(0)0,y ='(0)0,x ='(0) 1.049355751y =-(2)求解隐式微分方程组''''''''''''2235x y x y x y x y xy y ⎧+=⎨++-=⎩提示：运用符号盘算函数solve 求'''',x y ,然后运用求解微分方程的办法四．偏微分方程解法Matlab 供给了两种办法解决PDE 问题,一是运用pdepe 函数,它可以求解一般的PDEs,具有较大的通用性,但只支撑敕令情势挪用;二是运用PDE 对象箱,可以求解特别PDE 问题,PDEtoll 有较大的局限性,比如只能求解二阶PDE 问题,并且不能解决片微分方程组,但是它供给了GUI 界面,从庞杂的编程中摆脱出来,同时还可以经由过程File —>Save As 直接生成M 代码.1.一般偏微分方程（组）的求解(1)Matlab 供给的pdepe 函数,可以直接求解一般偏微分方程（组）,它的挪用格局为：sol=pdepe(m,@pdefun,@pdeic,@pdebc,x,t)@pdefun 是PDE 的问题描写函数,它必须换成标准情势：(,,)[(,,,)](,,,)m m u u u uc x t x x f x t u s x t u x t x x x -∂∂∂∂∂=+∂∂∂∂∂如许,PDE 就可以编写进口函数：[c,f,s]=pdefun(x,t,u,du),m,x,t 对应于式中相干参数,du 是u 的一阶导数,由给定的输入变量可表示出c,f,s 这三个函数.@pdebc 是PDE 的边界前提描写函数,它必须化为情势：(,,)(,,).*(,,,)0up x t u q x t u f x t u x ∂==∂于是边值前提可以编写函数描写为：[pa,qa,pb,qb]=pdebc(x,t,u,du),个中a 表示下边界,b 表示上边界.@pdeic 是PDE 的初值前提,必须化为情势：00(,)u x t u =,故可以运用函数描写为：u0=pdeic(x)sol 是一个三维数组,sol(:,:,i)表示i u 的解,换句话说,k u 对应x(i)和t(j)时的解为s ol(i,j,k),经由过程sol,我们可以运用pdeval 函数直接盘算某个点的函数值.(2)实例解释求解偏微分2111222221220.024()0.17()u u F u u t x u u F u u t x ⎧∂∂=--⎪⎪∂∂⎨∂∂⎪=+-⎪∂∂⎩个中, 5.7311.46()x x F x e e -=-且知足初始前提12(,0)1,(,0)0u x u x ==及边界前提1(0,)0,u t x ∂=∂221(0,)0,(1,)1,(1,)0u u t u t t x ∂===∂解：(1)对比给出的偏微分方程和pdepe 函数求解的标准情势,原方程改写为111221220.024()1.*()10.17u u F u u x u F u u u t x x ∂⎡⎤⎢⎥--⎡⎤⎡⎤⎡⎤∂∂∂=+⎢⎥⎢⎥⎢⎥⎢⎥-∂∂∂⎣⎦⎣⎦⎣⎦⎢⎥⎢⎥∂⎣⎦ 可见1121220.024()10,,,()10.17u F u u x m c f s F u u u x ∂⎡⎤⎢⎥--⎡⎤⎡⎤∂====⎢⎥⎢⎥⎢⎥-∂⎣⎦⎣⎦⎢⎥⎢⎥∂⎣⎦%目的PDE 函数function [c,f,s]=pdefun(x,t,u,du)c=[1;1];f=[0.024*du(1);0.17*du(2)];temp=u(1)-u(2);s=[-1;1].*(exp(5.73*temp)-exp(-11.46*temp))end(2)边界前提改写为：下边界2010.*00f u ⎡⎤⎡⎤⎡⎤+=⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦上边界1110.*000u f -⎡⎤⎡⎤⎡⎤+=⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦ %边界前提函数function [pa,qa,pb,qb]=pdebc(xa,ua,xb,ub,t)pa=[0;ua(2)];qa=[1;0];pb=[ub(1)-1;0];qb=[0;1];end(3)初值前提改写为：1210u u ⎡⎤⎡⎤=⎢⎥⎢⎥⎣⎦⎣⎦%初值前提函数function u0=pdeic(x)u0=[1;0];end(4)编写主调函数clcx=0:0.05:1;t=0:0.05:2;m=0;sol=pdepe(m,@pdefun,@pdeic,@pdebc,x,t);subplot(2,1,1)surf(x,t,sol(:,:,1))subplot(2,1,2)surf(x,t,sol(:,:,2))演习与思虑： This example illustrates the straightforward formulation, computation, and plotting of the solution of a single PDE.2()u u t x x π∂∂∂=∂∂∂This equation holds on an interval 01x ≤≤ for times 0t ≥. The PDE satisfies the initial condition (,0)sin u x x π= and boundary conditions(0,)0;(1,)0t u u t e t x π-∂=+=∂2.PDEtool 求解偏微分方程(1)PDEtool （GUI ）求解偏微分方程的一般步骤在Matlab 敕令窗口输入pdetool,回车,PDE 对象箱的图形用户界面(GUI)体系就启动了.从界说一个偏微分方程问题到完成解偏微分方程的定解,全部进程大致可以分为六个阶段Step 1 “Draw 模式”绘制平面有界区域Ω,经由过程公式把Matlab 体系供给的实体模子：矩形.圆.椭圆和多边形,组合起来,生成须要的平面区域.Step 2 “Boundary 模式”界说边界,声明不同边界段的边界前提.Step3 “PDE 模式”界说偏微分方程,肯定方程类型和方程系数c,a,f,d,依据具体情况,还可以在不同子区域声明不同系数.Step 4 “Mesh 模式”网格化区域Ω,可以掌握主动生成网格的参数,对生成的网格进行多次细化,使网格朋分更细更合理.Step 5 “Solve 模式”解偏微分方程,对于椭圆型方程可以激活并掌握非线性自顺应解题器来处理非线性方程;对于抛物线型方程和双曲型方程,设置初始边界前提后可以求出给准时刻t 的解;对于特点值问题,可以求出给定区间上的特点值.求解完成后,可以返回到Step 4,对网格进一步细化,进行再次求解.Step 6 “View 模式”盘算成果的可视化,可以经由过程设置体系供给的对话框,显示所求的解的表面图.网格图.等高线图和箭头梯形图.对于抛物线型和双曲线型问题的解还可以进行为画演示.(2)实例解释用法求解一个正方形区域上的特点值问题：12|0u u u u λ∂Ω⎧-∆-=⎪⎨⎪=⎩正方形区域为：11,1 1.x x -≤≤-≤≤(1)运用PDE 对象箱打开GUI 求解方程(2)进入Draw 模式,绘制一个矩形,然后双击矩形,在弹出的对话框中设置Left=-1,Bottom=-1,Width=2,Height=2,确认并封闭对话框(3)进入Boundary 模式,边界前提采用Dirichlet 前提的默认值(4)进入PDE 模式,单击对象栏PDE 按钮,在弹出的对话框中方程类型选择Eigenmodes,参数设置c=1,a=-1/2,d=1,确认后封闭对话框(5)单击对象栏的 按钮,对正方形区域进行初始网格剖分,然后再对网格进一步细化剖分一次(6)点开solve菜单,单击Parameters选项,在弹出的对话框中设置特点值区域为[ -20,20](7)单击Plot菜单的Parameters项,在弹出的对话框中选中Color.Height(3-D plot)和show mesh项,然后单击Done确认(8)单击对象栏的“=”按钮,开端求解。

如何使用MATLAB求解微分方程学习资料MATLAB是一种强大的数值计算和科学编程平台，可以用于求解微分方程和微分方程组。

在使用MATLAB求解微分方程之前，需要掌握一些基础知识，包括MATLAB的基本语法和常用的求解微分方程的技术。

下面是一些学习资料和步骤，帮助您使用MATLAB求解微分方程。

1.学习MATLAB基本语法和操作：首先，您需要学习MATLAB的基本语法和常用操作。

您可以参考MATLAB的官方文档、教程和手册，以及MATLAB的在线资源和视频教程。

这些资源可以帮助您掌握MATLAB的基本操作，建立良好的编程习惯。

2.学习求解微分方程的方法：在使用MATLAB求解微分方程之前，您需要了解一些常用的求解微分方程的方法，例如数值方法和解析方法。

数值方法包括欧拉法、龙格-库塔法、四阶龙格-库塔法等；解析方法包括分离变量法、线性微分方程的常系数齐次法和非齐次法等。

您可以参考微积分的教科书、在线资源和视频教程，掌握这些方法。

3. 使用MATLAB求解一阶微分方程：一阶微分方程是最简单的微分方程形式。

您可以首先尝试使用MATLAB求解一阶微分方程。

MATLAB提供了几个函数来求解一阶微分方程，例如ode45、ode23、ode113等。

您可以使用这些函数来解决特定的一阶微分方程，并观察结果。

可以使用plot 函数绘制微分方程的解，以获得更直观的理解。

4.使用MATLAB求解高阶微分方程：一旦您熟悉了使用MATLAB求解一阶微分方程的方法，您可以尝试使用同样的方法来求解高阶微分方程。

在求解高阶微分方程时，您需要将其转化为一组一阶微分方程。

例如，对于二阶线性微分方程，您可以引入一个新的变量来表示未知函数的导数，然后将其转化为一组一阶微分方程。

然后，您可以使用相同的求解函数来求解这组一阶微分方程。

5. 使用MATLAB求解微分方程组：对于多元微分方程组，MATLAB提供了更多的函数来求解。

例如，ode45s可以用于求解刚体动力学方程，ode23t可以用于求解刚体动力学方程。

第二讲 导数与微分方程一、实验内容1、实际引例（牛顿冷却模型）警察上午9点钟发现一被谋杀者，并测得尸体温度为32.4℃，一小时以后，尸体的温度变为31.7℃，尸体所在房间的温度是20℃。

如果人的正常体温为36.5℃，并知道热物体冷却速度与自身温度与外界温度之差成正比，试推断谋杀时间。

解设T 为尸体温度，t 为时间（十进制，如10.5代表10点30分）,温度的变化率（dt dT ）比例常数为K解题关键热物体冷却速度与自身温度与外界温度之差成正比。

即 dtdT =k*(T-20). 该方程就为微分方程，那么如何求解喃？以下为MATLAB 求解过程：T=dsolve('DT=k*(T-20)','T(9)=32.4') %T （9）表示在9点时的温度为32.4℃ ans =(62*exp(k*t))/(5*exp(9*k)) + 20由于10点室温为31.7℃，带入后可求K,命令如下：先输入 solve(''),然后把上面得到的表达式复制进去(62*exp(k*t))/(5*exp(9*k)) + 20k=solve('(62*exp(k*10))/(5*exp(9*k)) + 20=31.7') （带入t=10,T=31.7) k=-0.05810763080728074591965065204834T=(62*exp(k*t))/(5*exp(9*k)) + 20T=20.91915298056906312642973833139/exp(0.05810763080728074591965065204834*t)+20 T=vpa(T,6) (保留几位数)T =20.9192/exp(0.0581076*t) + 20.0t=solve('20.9192/exp(0.0581076*t) + 20.0=36.5')ans =4.08392390395064778268825821864340.083923*60=5.0354即谋杀时间 凌晨4点5分左右二、引例子所涉及的知识点1、首先要列出微分方程，简单的来说出现了变化率（dtdT ）。

第四讲 Matlab 求解微分方程（组）理论介绍：Matlab 求解微分方程（组）命令 求解实例：Matlab 求解微分方程（组）实例实际应用问题通过数学建模所归纳得到的方程，绝大多数都是微分方程，真正能得到代数方程的机会很少.另一方面，能够求解的微分方程也是十分有限的，特别是高阶方程和偏微分方程（组）.这就要求我们必须研究微分方程（组）的解法：解析解法和数值解法. 一．相关函数、命令及简介1.在Matlab 中，用大写字母D 表示导数，Dy 表示y 关于自变量的一阶导数，D2y 表示y 关于自变量的二阶导数，依此类推.函数dsolve 用来解决常微分方程（组）的求解问题，调用格式为：X=dsolve(‘eqn1’,’eqn2’,…)函数dsolve 用来解符号常微分方程、方程组，如果没有初始条件，则求出通解，如果有初始条件，则求出特解.注意，系统缺省的自变量为t2.函数dsolve 求解的是常微分方程的精确解法，也称为常微分方程的符号解.但是，有大量的常微分方程虽然从理论上讲，其解是存在的，但我们却无法求出其解析解，此时，我们需要寻求方程的数值解，在求常微分方程数值解方面，MATLAB 具有丰富的函数，我们将其统称为solver ，其一般格式为：[T,Y]=solver(odefun,tspan,y0)说明：(1)solver 为命令ode45、ode23、ode113、ode15s 、ode23s 、ode23t 、ode23tb 、ode15i 之一.(2)odefun 是显示微分方程'(,)y f t y =在积分区间tspan 0[,]f t t =上从0t 到f t 用初始条件0y 求解.(3)如果要获得微分方程问题在其他指定时间点012,,,,f t t t t 上的解，则令tspan 012[,,,]f t t t t =(要求是单调的).(4)因为没有一种算法可以有效的解决所有的ODE 问题，为此，Matlab 提供了多种求解器solver ，对于不同的ODE 问题，采用不同的solver.表1 Matlab中文本文件读写函数说明：ode23、ode45是极其常用的用来求解非刚性的标准形式的一阶微分方程（组）的初值问题的解的Matlab常用程序，其中：ode23采用龙格-库塔2阶算法，用3阶公式作误差估计来调节步长，具有低等的精度.ode45则采用龙格-库塔4阶算法，用5阶公式作误差估计来调节步长，具有中等的精度.3．在matlab命令窗口、程序或函数中创建局部函数时，可用内联函数inline，inline函数形式相当于编写M函数文件，但不需编写M-文件就可以描述出某种数学关系.调用inline函数，只能由一个matlab表达式组成，并且只能返回一个变量，不允许[u,v]这种向量形式.因而，任何要求逻辑运算或乘法运算以求得最终结果的场合，都不能应用inline函数，inline函数的一般形式为：FunctionName=inline(‘函数内容’, ‘所有自变量列表’)例如：（求解F(x)=x^2*cos(a*x)-b ,a,b是标量；x是向量）在命令窗口输入:Fofx=inline(‘x .^2*cos(a*x)-b’ , ‘x’,’a’,’b’); g= Fofx([pi/3 pi/3.5],4,1) 系统输出为：g=-1.5483 -1.7259注意：由于使用内联对象函数inline 不需要另外建立m 文件，所有使用比较方便，另外在使用ode45函数的时候，定义函数往往需要编辑一个m 文件来单独定义，这样不便于管理文件，这里可以使用inline 来定义函数. 二．实例介绍1.几个可以直接用Matlab 求微分方程精确解的实例 例1 求解微分方程2'2x y xy xe -+=程序：syms x y; y=dsolve(‘Dy+2*x*y=x*exp(-x^2)’,’x ’)例 2 求微分方程'0x xy y e +-=在初始条件(1)2y e =下的特解并画出解函数的图形.程序：syms x y; y=dsolve(‘x*Dy+y-exp(1)=0’,’y(1)=2*exp(1)’,’x ’);ezplot(y)例 3 求解微分方程组530tdx x y e dtdy x y dt⎧++=⎪⎪⎨⎪--=⎪⎩在初始条件00|1,|0t t x y ====下的特解并画出解函数的图形.程序：syms x y t[x,y]=dsolve('Dx+5*x+y=exp(t)','Dy-x-3*y=0','x(0)=1','y(0)=0','t') simple(x); simple(y)ezplot(x,y,[0,1.3]);axis auto2.用ode23、ode45等求解非刚性标准形式的一阶微分方程（组）的初值问题的数值解（近似解）例 4 求解微分方程初值问题2222(0)1dy y x xdx y ⎧=-++⎪⎨⎪=⎩的数值解，求解范围为区间[0,0.5].程序：fun=inline('-2*y+2*x^2+2*x','x','y');[x,y]=ode23(fun,[0,0.5],1); plot(x,y,'o-')例 5 求解微分方程22'2(1)0,(0)1,(0)0d y dyy y y y dt dtμ--+===的解，并画出解的图形.分析：这是一个二阶非线性方程，我们可以通过变换，将二阶方程化为一阶方程组求解.令12,,7dyx y x dtμ===，则 121221212,(0)17(1),(0)0dx x x dtdx x x x x dt⎧==⎪⎪⎨⎪=--=⎪⎩ 编写M-文件vdp.m function fy=vdp(t,x)fy=[x(2);7*(1-x(1)^2)*x(2)-x(1)]; end在Matlab 命令窗口编写程序 y0=[1;0][t,x]=ode45(@vdp,[0,40],y0);或[t,x]=ode45('vdp',[0,40],y0); y=x(:,1);dy=x(:,2); plot(t,y,t,dy)练习与思考：M-文件vdp.m 改写成inline 函数程序？ 3.用Euler 折线法求解Euler 折线法求解的基本思想是将微分方程初值问题00(,)()dyf x y dxy x y ⎧=⎪⎨⎪=⎩ 化成一个代数(差分)方程，主要步骤是用差商()()y x h y x h +-替代微商dydx，于是00()()(,())()k k k k y x h y x f x y x h y y x +-⎧=⎪⎨⎪=⎩记1,(),k k k k x x h y y x +=+=从而1(),k k y y x h +=+于是0011(),,0,1,2,,1(,).k k k k k k y y x x x h k n y y hf x y ++=⎧⎪=+=-⎨⎪=+⎩例 6 用Euler 折线法求解微分方程初值问题22(0)1dyx y dxy y ⎧=+⎪⎨⎪=⎩的数值解（步长h 取0.4），求解范围为区间[0,2].分析：本问题的差分方程为00110,1,0.4,0,1,2,,1(,).k k k k k k x y h x x h k n y y hf x y ++===⎧⎪=+=-⎨⎪=+⎩程序：>> clear >> f=sym('y+2*x/y^2'); >> a=0; >> b=2; >> h=0.4; >> n=(b-a)/h+1; >> x=0; >> y=1;>> szj=[x,y];%数值解 >> for i=1:n-1y=y+h*subs(f,{'x','y'},{x,y});%subs ，替换函数 x=x+h; szj=[szj;x,y]; end >>szj>> plot(szj(:,1),szj(:,2))说明：替换函数subs 例如：输入subs(a+b,a,4) 意思就是把a 用4替换掉，返回 4+b ，也可以替换多个变量，例如：subs(cos(a)+sin(b),{a,b},[sym('alpha'),2])分别用字符alpha 替换a 和2替换b ，返回 cos(alpha)+sin(2)特别说明：本问题可进一步利用四阶Runge-Kutta 法求解，Euler 折线法实际上就是一阶Runge-Kutta 法，Runge-Kutta 法的迭代公式为001112341213243(),,(22),6(,),0,1,2,,1(,),22(,),22(,).k k k k k k k k k k k k y y x x x h h y y L L L L L f x y k n h h L f x y L h h L f x y L L f x h y hL ++=⎧⎪=+⎪⎪=++++⎪⎪=⎪=-⎨⎪=++⎪⎪⎪=++⎪⎪=++⎩相应的Matlab 程序为：>> clear >> f=sym('y+2*x/y^2'); >> a=0; >> b=2; >> h=0.4; >> n=(b-a)/h+1; >> x=0; >> y=1;>> szj=[x,y];%数值解 >> for i=1:n-1l1=subs(f, {'x','y'},{x,y});替换函数 l2=subs(f, {'x','y'},{x+h/2,y+l1*h/2}); l3=subs(f, {'x','y'},{x+h/2,y+l2*h/2}); l4=subs(f, {'x','y'},{x+h,y+l3*h}); y=y+h*(l1+2*l2+2*l3+l4)/6; x=x+h; szj=[szj;x,y]; end>>szj>> plot(szj(:,1),szj(:,2))练习与思考：(1)ode45求解问题并比较差异. (2)利用Matlab 求微分方程(4)(3)''20y y y -+=的解.(3)求解微分方程''2',2(1)0,030,(0)1,(0)0y y y y x y y --+=≤≤==的特解. (4)利用Matlab 求微分方程初值问题2''''00(1)2,|1,|3x x x y xy y y ==+===的解. 提醒：尽可能多的考虑解法 三．微分方程转换为一阶显式微分方程组Matlab 微分方程解算器只能求解标准形式的一阶显式微分方程（组）问题，因此在使用ODE 解算器之前，我们需要做的第一步，也是最重要的一步就是借助状态变量将微分方程（组）化成Matlab 可接受的标准形式.当然，如果ODEs 由一个或多个高阶微分方程给出，则我们应先将它变换成一阶显式常微分方程组.下面我们以两个高阶微分方程组构成的ODEs 为例介绍如何将它变换成一个一阶显式微分方程组.Step 1 将微分方程的最高阶变量移到等式左边，其它移到右边，并按阶次从低到高排列.形式为：()'''(1)'''(1)()'''(1)'''(1)(,,,,,,,,,,)(,,,,,,,,,,)m m n n m n x f t x x x x y y y y y g t x x x x y y y y ----⎧=⎨=⎩Step 2 为每一阶微分式选择状态变量，最高阶除外'''(1)123'''(1)123,,,,,,,,,m m n m m m m n x x x x x x x x x y x y x y x y--++++========注意：ODEs 中所有是因变量的最高阶次之和就是需要的状态变量的个数，最高阶的微分式不需要给它状态变量.Step 3 根据选用的状态变量，写出所有状态变量的一阶微分表达式''''122334123''12123,,,,(,,,,,),,(,,,,,)m m n m m m nm n x x x x x x x f t x x x x xx xg t x x x x +++++======练习与思考：(1)求解微分方程组**'''3312*'''3312()()22x x x y x r r y y y x y r r μμμμμμ⎧+-=+--⎪⎪⎨⎪=+--⎪⎩其中2r =1r =*1,μμ=-1/82.45,μ=(0) 1.2,x =(0)0,y ='(0)0,x ='(0) 1.049355751y =-(2)求解隐式微分方程组''''''''''''2235x y x y x y x y xy y ⎧+=⎨++-=⎩ 提示：使用符号计算函数solve 求'''',x y ，然后利用求解微分方程的方法 四．偏微分方程解法Matlab 提供了两种方法解决PDE 问题，一是使用pdepe 函数，它可以求解一般的PDEs,具有较大的通用性，但只支持命令形式调用；二是使用PDE 工具箱，可以求解特殊PDE 问题，PDEtoll 有较大的局限性，比如只能求解二阶PDE 问题，并且不能解决片微分方程组，但是它提供了GUI 界面，从复杂的编程中解脱出来，同时还可以通过File —>Save As 直接生成M 代码.1.一般偏微分方程（组）的求解(1)Matlab 提供的pdepe 函数，可以直接求解一般偏微分方程（组），它的调用格式为：sol=pdepe(m,@pdefun,@pdeic,@pdebc,x,t)@pdefun 是PDE 的问题描述函数，它必须换成标准形式：(,,)[(,,,)](,,,)m m u u u uc x t x x f x t u s x t u x t x x x-∂∂∂∂∂=+∂∂∂∂∂ 这样，PDE 就可以编写入口函数：[c,f,s]=pdefun(x,t,u,du)，m,x,t 对应于式中相关参数，du 是u 的一阶导数，由给定的输入变量可表示出c,f,s 这三个函数.@pdebc 是PDE 的边界条件描述函数，它必须化为形式：(,,)(,,).*(,,,)0up x t u q x t u f x t u x∂==∂ 于是边值条件可以编写函数描述为：[pa,qa,pb,qb]=pdebc(x,t,u,du),其中a 表示下边界，b 表示上边界.@pdeic 是PDE 的初值条件，必须化为形式：00(,)u x t u =，故可以使用函数描述为：u0=pdeic(x)sol 是一个三维数组，sol(:,:,i)表示i u 的解，换句话说，k u 对应x(i)和t(j)时的解为sol(i,j,k)，通过sol ，我们可以使用pdeval 函数直接计算某个点的函数值.(2)实例说明 求解偏微分2111222221220.024()0.17()u u F u u t xu u F u u tx ⎧∂∂=--⎪⎪∂∂⎨∂∂⎪=+-⎪∂∂⎩ 其中， 5.7311.46()x x F x e e -=-且满足初始条件12(,0)1,(,0)0u x u x ==及边界条件1(0,)0,u t x ∂=∂221(0,)0,(1,)1,(1,)0uu t u t t x∂===∂ 解：(1)对照给出的偏微分方程和pdepe 函数求解的标准形式，原方程改写为111221220.024()1.*()10.17u u F u u x u F u u u t x x ∂⎡⎤⎢⎥--⎡⎤⎡⎤⎡⎤∂∂∂=+⎢⎥⎢⎥⎢⎥⎢⎥-∂∂∂⎣⎦⎣⎦⎣⎦⎢⎥⎢⎥∂⎣⎦可见1121220.024()10,,,()10.17u F u u x m c f s F u u u x ∂⎡⎤⎢⎥--⎡⎤⎡⎤∂====⎢⎥⎢⎥⎢⎥-∂⎣⎦⎣⎦⎢⎥⎢⎥∂⎣⎦ %目标PDE 函数function [c,f,s]=pdefun(x,t,u,du) c=[1;1];f=[0.024*du(1);0.17*du(2)]; temp=u(1)-u(2);s=[-1;1].*(exp(5.73*temp)-exp(-11.46*temp)) end(2)边界条件改写为：下边界2010.*00f u ⎡⎤⎡⎤⎡⎤+=⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦上边界1110.*000u f -⎡⎤⎡⎤⎡⎤+=⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦%边界条件函数function [pa,qa,pb,qb]=pdebc(xa,ua,xb,ub,t) pa=[0;ua(2)]; qa=[1;0]; pb=[ub(1)-1;0]; qb=[0;1]; end(3)初值条件改写为：1210u u ⎡⎤⎡⎤=⎢⎥⎢⎥⎣⎦⎣⎦%初值条件函数 function u0=pdeic(x) u0=[1;0]; end(4)编写主调函数 clc x=0:0.05:1; t=0:0.05:2; m=0;sol=pdepe(m,@pdefun,@pdeic,@pdebc,x,t); subplot(2,1,1) surf(x,t,sol(:,:,1)) subplot(2,1,2) surf(x,t,sol(:,:,2))练习与思考： This example illustrates the straightforward formulation, computation, and plotting of the solution of a single PDE.2()u u t x xπ∂∂∂=∂∂∂ This equation holds on an interval 01x ≤≤ for times 0t ≥. The PDE satisfies the initial condition (,0)sin u x x π= and boundary conditions(0,)0;(1,)0t uu t e t xπ-∂=+=∂ 2.PDEtool 求解偏微分方程(1)PDEtool （GUI ）求解偏微分方程的一般步骤在Matlab 命令窗口输入pdetool ，回车，PDE 工具箱的图形用户界面(GUI)系统就启动了.从定义一个偏微分方程问题到完成解偏微分方程的定解，整个过程大致可以分为六个阶段Step 1 “Draw 模式”绘制平面有界区域Ω，通过公式把Matlab 系统提供的实体模型：矩形、圆、椭圆和多边形，组合起来，生成需要的平面区域.Step 2 “Boundary 模式”定义边界，声明不同边界段的边界条件.Step 3 “PDE 模式”定义偏微分方程，确定方程类型和方程系数c,a,f,d ，根据具体情况，还可以在不同子区域声明不同系数.Step 4 “Mesh 模式”网格化区域Ω，可以控制自动生成网格的参数，对生成的网格进行多次细化，使网格分割更细更合理.Step 5 “Solve 模式”解偏微分方程，对于椭圆型方程可以激活并控制非线性自适应解题器来处理非线性方程；对于抛物线型方程和双曲型方程，设置初始边界条件后可以求出给定时刻t 的解；对于特征值问题，可以求出给定区间上的特征值.求解完成后，可以返回到Step 4，对网格进一步细化，进行再次求解.Step 6 “View 模式”计算结果的可视化，可以通过设置系统提供的对话框，显示所求的解的表面图、网格图、等高线图和箭头梯形图.对于抛物线型和双曲线型问题的解还可以进行动画演示.(2)实例说明用法求解一个正方形区域上的特征值问题：12|0u u u u λ∂Ω⎧-∆-=⎪⎨⎪=⎩ 正方形区域为：11,1 1.x x -≤≤-≤≤(1)使用PDE 工具箱打开GUI 求解方程(2)进入Draw 模式，绘制一个矩形，然后双击矩形，在弹出的对话框中设置Left=-1,Bottom=-1,Width=2,Height=2,确认并关闭对话框(3)进入Boundary 模式，边界条件采用Dirichlet 条件的默认值(4)进入PDE 模式，单击工具栏PDE 按钮，在弹出的对话框中方程类型选择Eigenmodes,参数设置c=1,a=-1/2,d=1,确认后关闭对话框(5)单击工具栏的 按钮，对正方形区域进行初始网格剖分，然后再对网格进一步细化剖分一次(6)点开solve菜单，单击Parameters选项，在弹出的对话框中设置特征值区域为[-20,20](7)单击Plot菜单的Parameters项，在弹出的对话框中选中Color、Height(3-D plot)和show mesh项，然后单击Done确认(8)单击工具栏的“=”按钮，开始求解。

第五讲 Matlab求解微分方程教学目的：学会用MATLAB求简单微分方程的解析解、数值解，加深对微分方程概念和应用的理解；针对一些具体的问题，如追击问题，掌握利用软件求解微分方程的过程；了解微分方程模型解决问题思维方法及技巧；体会微分方程建摸的艺术性．教学重点：利用机理分析建模微分方程模型，掌握追击问题的建模方法，掌握利用MATLAB求解数值解．教学难点：利用机理分析建模微分方程模型，通过举例，结合图形以及与恰当的假设突破教学难点．1微分方程相关函数（命令）及简介因为没有一种算法可以有效地解决所有的ODE 问题，为此，Matlab 提供了多种求解器Solver，对于不同的ODE 问题，采用不同的Solver．阶常微分方程(组)的初值问题的解的 Matlab 的常用程序，其中：ode23 采用龙格-库塔2 阶算法，用3 阶公式作误差估计来调节步长，具有低等的精度．ode45 则采用龙格-库塔4 阶算法，用5 阶公式作误差估计来调节步长，具有中等的精度．2 求解微分方程的一些例子2.1 几个可以直接用 Matlab 求微分方程精确解的例子：例1：求解微分方程22x xe xy dxdy-=+，并加以验证． 求解本问题的Matlab 程序为：syms x y %line1 y=dsolve('Dy+2*x*y=x*exp(-x^2)','x') %line2 diff(y,x)+2 *x*y-x*exp(-x^2) %line3 simplify(diff(y,x)+2*x*y-x*exp(-x^2)) %line4说明：(1) 行line1是用命令定义x,y 为符号变量．这里可以不写，但为确保正确性，建议写上；(2) 行line2是用命令求出的微分方程的解：1/2*exp(-x^2)*x^2+exp(-x^2)*C1(3) 行line3使用所求得的解．这里是将解代入原微分方程，结果应该为0，但这里给出：-x^3*exp(-x^2)-2*x*exp(-x^2)*C1+2*x*(1/2*exp(-x^2)*x^2+exp(-x^2)*C1)(4) 行line4 用 simplify() 函数对上式进行化简，结果为 0， 表明)(x y y =的确是微分方程的解．例2：求微分方程0'=-+x e y xy 在初始条件e y 2)1(=下的特解，并画出解函数的图形．求解本问题的 Matlab 程序为： syms x yy=dsolve('x*Dy+y-exp(x)=0','y(1)=2*exp(1)','x') ezplot(y)微分方程的特解为：y=1/x*exp(x)+1/x* exp (1) (Matlab 格式)，即xe e y x+=，此函数的图形如图 1：图1 y 关于x 的函数图象2.2 用ode23、ode45等求解非刚性的标准形式的一阶常微分方程(组)的初值问题的数值解(近似解)．例３：求解微分方程初值问题⎪⎩⎪⎨⎧=++-=1)0(2222y xx y dx dy 的数值解，求解范围为区间[0, 0.5]．fun=inline('-2*y+2*x^2+2*x','x','y'); [x,y]=ode23(fun,[0,0.5],1); x'; y'; plot(x,y,'o-') >> x' ans =0.0000 0.0400 0.0900 0.1400 0.1900 0.2400 0.2900 0.3400 0.3900 0.4400 0.4900 0.5000>> y' ans =1.0000 0.9247 0.8434 0.7754 0.7199 0.6764 0.6440 0.6222 0.6105 0.6084 0.6154 0.6179 图形结果为图 2．图2 y 关于x 的函数图像3 常微分在实际中的应用 3.1 导弹追踪问题1、符号说明，w ，乙舰的速率恒为v 0；设时刻t 乙舰的坐标为((),())X t Y t ，导弹的坐标为((),())x t y t ；当零时刻，((0),(0))(1,0)X Y =，((0),(0))(0,0)x y =．建立微分方程模型．202,01(0)0,'(0)0v d yk x dx w y y ⎧⎪⎪= =<<⎨⎪⎪==⎩由微分方程模型解得11(1)(1)11(),12(1)2(1)2(1)2(1)k k x x y x k k k k k +-+--=-+-≠+-+-++代入题设的数据1/5k =，得到导弹的运行轨迹为4655555(1)(1)81224y x x =--+-+当1=x 时245=y ，即当乙舰航行到点)245,1(处时被导弹击中. 被击中时间为:00245v v y t ==. 若v 0=1, 则在t =0.21处被击中. 利用MALAB 作图如图３．clear, x=0:0.01:1; y=-5*(1-x).^(4/5)/8+5*(1-x).^(6/5)/12+5/24; plot(x,y,'*')图３导弹运行轨迹（解析法） 图４两种方法对比的导弹运行轨迹2、数值方法求解．设导弹速率恒为w ，则得到参数方程为⎪⎪⎩⎪⎪⎨⎧--+-=--+-=)()()()()()(2222y Y y Y x X wdt dy x X y Y x X w dt dx因乙舰以速度0v 沿直线1x =运动，设01v =，5,1,w X Y t ===，因此导弹运动轨迹的参数方程为：⎪⎪⎪⎩⎪⎪⎪⎨⎧==--+-=--+-=0)0(,0)0()()()1(5)1()()1(52222y x y t y t x dtdyx y t x dt dxMATLAB 求解数值解程序如下，结果见图４ t0=0,tf=0.21;[t,y]=ode45('eq2',[t0 tf],[0 0]); X=1;Y=0:0.001:0.21;plot(X,Y,'-') plot(y(:,1),y(:,2),'*'),hold onx=0:0.01:1; y=-5*(1-x).^(4/5)/8+5*(1-x).^(6/5)/12+5/24; plot(x,y,'r')很明显数值计算的方法比较简单而且适用． 3.2 蚂蚁追击问题（1）建立平面直角坐标系．以正多边形的中心为原点, 设正多边形的一个顶点为起始点, 连接此点和原点作x 轴. 根据x 轴作出相应的y 轴; 选取足够小的t ∆进行采样．（2）在每一时刻t ，计算每只蚂蚁在下一时刻t t +∆时的坐标．不妨设甲追逐对象是乙，在时间t 时，甲的坐标为A 11(,)x y ，乙的坐标为B 22(,)x y ．甲在t t +∆时在'A 点(如图1), 其坐标为11(cos ,sin )x v t y v t θθ+∆+∆，其中2221212121cos ,sin ,()()x x y yd x x y y d dθθ--===-+-. 同理，依次计算下一只蚂蚁在t t +∆时的坐标．通过间隔t ∆进行采样，得到新一轮各个蚂蚁在一个新的正多边形位置坐标．（4）重复2）步，直到d 充分小为止．（5）连接每只蚂蚁在各时刻的位置，就得到所求的轨迹．用MALAB 求解并作图，函数zhuJi(x,y)在附录一定义，如图６ t=[1:8]; s=7*exp(t.*2*pi/length(t)*i); x=real(s); y=imag(s);zhuJi(x,y)图6当蚂蚁为7只时的图形习题1. 求微分方程0sin 2')1(2=-+-x xy y x 的通解．2. 求微分方程x e y y y x sin 5'2''=+-的通解．3. 求微分方程组'A 11,(,)A x y 22,(,)B x y dθ图1 在采样时间内，相连蚂蚁追击⎪⎪⎩⎪⎪⎨⎧=-+=++00y x dtdy y x dtdx在初始条件0|,1|00====t t y x 下的特解，并画出解函数()y f x =的图形．4. 分别用 ode23、ode45 求上述第 3 题中的微分方程初值问题的数值解(近似解)，求解区间为[0,2]t ∈．利用画图来比较两种求解器之间的差异．4 参考文献:[1] Mastering Matlab 6, D. Hanselman, B. Littlefield, 清华大学出版,2002[2] 赵静等编．数学建模与数学实验（第3版）．北京：高等教育出版社．2008. [3] 姜启源编. 数学模型（第二版）．北京：高等教育出版社.1993.[4] 石勇国. 蚂蚁追击问题与等角螺线. 宜宾学院学报. 2008，(6): 23-25. [5] 张伟年，杜正东，徐冰．常微分方程．北京：高等教育出版社．2006．5 附录附录一:zhuji(x,y)的M 文件 function zhuji(x,y) clf v=1; dt=0.05;x(length(x)+1)=x(1);y(length(y)+1)=y(1); plot(x,y,'*') holdfor it=1:1000for i=1:length(x)-1d=sqrt((x(i)-x(i+1))^2+(y(i)-y(i+1))^2); x(i)=x(i)+v*dt*(x(i+1)-x(i))/d; y(i)=y(i)+v*dt*(y(i+1)-y(i))/d; endplot(x,y,'.') hold onx(length(x))=x(1); y(length(y))=y(1); end。