1、 1 第四讲 Matlab 求解微分方程（组） 理论介绍： Matlab 求解微分方程（组）命令 求解实例： Matlab 求解微分方程（组）实例 实际应用问题通过数学建模所归纳得到的方程，绝大多数都是微分方程，真正能得到代数方程的机会很少 .另一方面，能够求解的微分方程也是十分有限的，特别是高阶方程和偏微分方程（组） .这就要求我们必须研究微分方程（组）的解法：解析解法和数值解法 . 一 相关函数、命令及简介 1.在 Matlab 中，用大写字母 D 表示导数， Dy 表示 y 关于自变量的一阶导数，D2y 表示 y 关于自变量的二阶导数，依此类推 .函数 dsolve 用来解决常微分方程（
2、组）的求解问题，调用格式为： X=dsolve(eqn1,eqn2,) 函数 dsolve 用来解符号常微分方程、方程组，如果没有初始条件，则求出通解，如果有初始条件，则求出特解 . 注意 ， 系统缺省的自变量为 t 2.函数 dsolve 求解的是 常微分方程的精确解法，也称为常微分方程的符号解 .但是，有大量的常微分方程虽然从理论上讲，其解是存在的，但我们却无法求出其解析解，此时，我们需要寻求方程的数值解，在求常微分方程数值解方面，MATLAB 具有丰富的函数，我们将其统称为 solver，其一般格式为： T,Y=solver(odefun,tspan,y0) 说明 ： (1)solver
3、 为命令 ode45、 ode23、 ode113、 ode15s、 ode23s、 ode23t、 ode23tb、ode15i 之一 . (2)odefun 是 显示微分方程 ( , )y f t y 在积分区间 tspan 0 , ftt 上从 0t 到 ft 用初 始条件 0y 求解 . (3)如果要获得微分方程问题在其他指定时间点 0 1 2, , , , ft t t t 上的解，则令tspan 0 1 2 , , , ft t t t (要求是单调的 ). (4)因为没有一种算法可以有效的解决所有的 ODE 问题，为此， Matlab 提供了多种求解器 solver，对于不同的
4、ODE 问题，采用不同的 solver. 2 表 1 Matlab 中文本文件读写函数 求解器 ODE 类型 特点 说明 ode45 非刚性 单步算法 ： 4、 5 阶 Runge-Kutta方程；累计截断误差 3()x 大部分场合的首选算法 ode23 非刚性 单步算法： 2、 3 阶 Runge-Kutta方程；累计截断误差 3()x 使用于精度较低的情形 ode113 非刚性 多步法： Adams 算法；高低精度可达 3610 10 计算时间比 ode45短 ode23t 适度刚性 采用梯形算法 适度刚性情形 ode15s 刚性 多步法： Gears 反向数值微分；精度中等 若 ode4
5、5 失效时，可尝试使用 ode23s 刚性 单步法： 2 阶 Rosebrock 算法；低精度 当精度较低时，计算时间比 ode15s 短 ode23tb 刚性 梯形算法；低精度 当精度较低时，计算时间比 ode15s 短 说明 ： ode23、 ode45 是极其常用的用来求解非刚性的标准形式的一阶微分方程（组）的初值问题的解的 Matlab 常用程序，其中： ode23 采用龙格 -库塔 2 阶算法，用 3 阶公式作误差估计来调节步长，具有低等的精度 . ode45 则采用龙格 -库塔 4 阶算法，用 5 阶 公式作误差估计来调节步长，具有中等的精度 . 3 在 matlab 命令窗口、程
6、序或函数中创建局部函数时，可用 内联函数 inline，inline 函数 形式相当于编写 M 函数 文件 ，但不需编写 M-文件就可以描述出某种数学关系 .调用 inline 函数，只能由一个 matlab 表达式组成，并且只能返回一个变量 ， 不允许 u,v这种 向量 形式 .因而，任何要求逻辑运算或乘法运算以求得最终结果的场合，都不能应用 inline 函数 ， inline 函数的一般形式为： FunctionName=inline(函数内容 , 所有自变量列表 ) 例 如：（求解 F(x)=x2*cos(a*x)-b ,a,b 是标量； x 是向量 ）在命令窗口输入 : 3 Fofx
7、=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 2 xy xy xe 程序： syms x y; y=dsolve(Dy+2*x*y=x*exp
8、(-x2),x) 例 2 求微分方程 0xxy y e 在初始条件 (1) 2ye 下的特解 并画出解函数的图形 . 程序： syms x y; y=dsolve(x*Dy+y-exp(1)=0,y(1)=2*exp(1),x);ezplot(y) 例 3 求解微分方程组530tdx x y edtdy xydt 在初始条件 00| 1, | 0ttxy下的特解并画出解函数的图形 . 程序： 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.
9、3);axis auto 2.用 ode23、 ode45 等求解非刚性标准形式 的一阶微分方程（组）的初值问题的数值解（近似解） 例 4 求解微分方程初值问题 22 2 2(0 ) 1dy y x xdxy 的数值解，求解范围为区间 0,0.5. 程序： fun=inline(-2*y+2*x2+2*x,x,y); 4 x,y=ode23(fun,0,0.5,1); plot(x,y,o-) 例 5 求解微分方程 2 22 ( 1 ) 0 , (0 ) 1 , (0 ) 0d y d yy y y yd t d t 的解，并画出解的图形 . 分析：这是一个二阶非线性方程，我们可以通过变换，将
10、二阶方程化为一阶方程组 求解 .令12, , 7dyx y x dt ，则 121221 2 1 2, ( 0 ) 17 (1 ) , ( 0 ) 0dx xxdtdx x x x xdt 编写 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 改写成
11、inline 函数程序？ 3.用 Euler 折线法求解 Euler 折线法求解的基本思想是将微分方程初值问题 00( , )()dy f x ydxy x y 化成一个代数 (差分 )方程，主要步骤是用差商 ( ) ( )y x h y xh 替代微商 dydx ，于是 00( ) ( ) ( , ( ) )()kk kky x h y x f x y xhy y x 5 记 1 , ( ),k k k kx x h y y x 从而 1 ( ),kky x h 于是 0011( ) , 0 , 1 , 2 , , 1( , ) .kkk k k ky y xx x h k ny y h f
12、 x y 例 6 用 Euler 折线法 求解微分方程初值问题 22(0) 1dy xydx yy 的数值解 （ 步长 h 取 0.4），求解范围为区间 0,2. 分析： 本 问题的差分方程为 00110 , 1 , 0 .4, 0 , 1 , 2 , , 1( , ) .kkk k k kx y hx x h k ny y h f x y 程序： clear f=sym(y+2*x/y2); a=0; b=2; h=0.4; n=(b-a)/h+1; x=0; y=1; szj=x,y;%数值解 for i=1:n-1 y=y+h*subs(f,x,y,x,y);%subs，替换函数 x=x
13、+h; szj=szj;x,y; end szj plot(szj(:,1),szj(:,2) 说明 ：替换函数 subs 例如 ： 输入 subs(a+b,a,4) 意思就是把 a 用 4 替换掉，6 返回 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 法的迭代公式为 0011 1
14、 2 3 41213243( ) ,( 2 2 ) ,6( , ) ,0 , 1 , 2 , , 1( , ) ,22( , ) ,22( , ) .kkkkkkkkkkkky y xx x hhy y L L L LL f x yknhhL f x y LhhL f x y LL f x h y hL 相应的 Matlab 程序为： clear f=sym(y+2*x/y2); a=0; b=2; h=0.4; n=(b-a)/h+1; x=0; y=1; szj=x,y;%数值解 for i=1:n-1 l1=subs(f, x,y,x,y);替换函数 l2=subs(f, x,y,x+h
15、/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 7 szj plot(szj(:,1),szj(:,2) 练习与思考 ： (1)ode45 求解问题并比较差异 . (2)利用 Matlab 求 微分方程 ( 4 ) (3) 20y y y 的解 . (3)求解微分方程 2 ,2 ( 1 ) 0 , 0 3 0 , ( 0 ) 1 , ( 0 ) 0y y y y x y y 的特解 . (4)利
16、用 Matlab 求 微分方程 初值问题 2 00(1 ) 2 , | 1 , | 3xxx y x y y y 的解 . 提醒 ：尽可能多的考虑解法 三 微分方程转换为一阶显式微分方程组 Matlab 微分方程解算器只能求解标准形式的一阶显式微分方程（组）问题，因此在使用 ODE 解算器之前，我们需要做的第一步，也是最重要的一步就是借助状态变量将微分方程（组）化成 Matlab 可接受的标准形式 .当然，如果 ODEs由一个或多个高阶微分方程给出，则我们应先将它变换成一阶显式常微分方程组 .下面我们以两个高阶微分方程组构成的 ODEs 为例介绍如何将它变换成一个一阶显式微分方程组 . Ste
17、p 1 将微分方程的最高阶变量移到等式左边，其它移到右边，并按阶次从低到高排列 .形式为： ( ) ( 1 ) ( 1 )( ) ( 1 ) ( 1 )( , , , , , , , , , , )( , , , , , , , , , , )m m nn m nx f t x x x x y y y yy g t x x x x y y y y Step 2 为每一阶微分式选择状态变量，最高阶除外 ( 1 )1 2 3 ( 1 )1 2 3, , , , , , , , mmnm m m m nx x x x x x x xx y x y x y x y 注意 ： ODEs 中所有是因变量的
18、最高阶次之和就是需要的状态变量的个数，最高阶的微分式不需要给它状态变量 . Step 3 根据选用的状态变量，写出所有状态变量的一阶微分表达式 1 2 2 3 3 4 1 2 31 2 1 2 3, , , , ( , , , , , ), , ( , , , , , )m m nm m m n m nx x x x x x x f t x x x xx x x g t x x x x 练习与思考 ： (1)求解微分方程组 8 * 3312* 3312( ) ( )22xxx y xrryyy x yrr 其中 * 2 22 ( ) ,r x y 221 ( ) ,r x y * 1, 1/8
19、2.45, (0) 1.2,x (0) 0,y (0) 0,x (0 ) 1 .0 4 9 3 5 5 7 5 1y (2)求解隐式微分方程组 2235x y x yx y x y xy y 提示 ：使用符号计算函数 solve 求 ,xy，然后利用求解微分方程的方法 四 偏 微分方程 解法 Matlab 提供了两种方法解决 PDE 问题，一是使用 pdepe 函数，它可以求解一般的 PDEs,具有较大的通用性，但只支持命令形式调用；二是使用 PDE 工具箱，可以求解特殊 PDE 问题， PDEtoll 有较大的局限性，比如只能求解二阶 PDE问题，并且不能解决片微分方程组，但是它提供了 GU
20、I 界面，从复杂的编程中解脱出来， 同时还可以通过 File Save As 直接生成 M 代码 . 1.一般偏微分方程（组）的求解 (1)Matlab 提 供的 pdepe 函数，可以直接求解一般偏微分方程（组），它的调用格式为： sol=pdepe(m,pdefun,pdeic,pdebc,x,t) pdefun 是 PDE 的问题描述函数，它必须换成标准形式： ( , , ) ( , , , ) ( , , , )mmu u u uc x t x x f x t u s x t ux t x x x 这样， PDE 就可以编写入口函数： c,f,s=pdefun(x,t,u,du)， m
21、,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 ，故可以使用函数9 描述为： u0=pdeic(x) sol 是一个三维数组， sol(:,:
22、,i)表示 iu 的解 ，换句话说， ku 对应 x(i)和 t(j)时的解为 sol(i,j,k)，通过 sol，我们可以使用 pdeval 函数直接计算某个点的函数值 . (2)实例说明 求解偏微分 2111222221220 .0 2 4 ( )0 .1 7 ( )uu F u utxuu F u utx 其中， 5 .7 3 1 1 .4 6() xxF x e e且满足初始条件 12( , 0 ) 1, ( , 0 ) 0u x u x及边界条件1 (0, ) 0,u tx 221(0 , ) 0 , (1 , ) 1 , (1 , ) 0uu t u t tx 解 ： (1)对照给
23、出的偏微分方程 和 pdepe 函数求解的标准形 式 ，原方程改写为 11 1 22 1 220 .0 2 4 ()1.* ()10 .1 7uu F u uxu F u uutxx 可见1121220 .0 2 4 ()10 , , , ()10 .1 7uF u uxm c f sF u uux %目标 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)边界条件改写为：
24、下边界20 10.*00fu 上边界 1 1 1 0.*0 0 0u f 10 %边界条件函数 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)初值条件改写为： 1210uu%初值条件函数 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,s
25、ol(:,:,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 ()uut x x This equation holds on an interval 01x for times 0t . The PDE satisfies the initial condition ( ,0) sinu x x and bound
26、ary conditions (0 , ) 0 ; (1 , ) 0t uu t e tx 2.PDEtool 求解偏微分方程 11 (1)PDEtool（ GUI）求 解偏微分方程的一般步骤 在 Matlab 命令窗口输入 pdetool，回车， PDE 工具箱的图形用户界面 (GUI)系统就启动了 .从定义一个偏微分方程问题到完成解偏微分方程的定解，整个过程大致可以分为六个阶段 Step 1 “ Draw 模式”绘制平面有界区域 ，通过公式把 Matlab 系统提供的实体模型：矩形、圆、椭圆和多边形，组合起来，生成需要的平面区域 . Step 2 “ Boundary 模式”定义边界，声明
27、不同边界段的边界条件 . Step 3 “ PDE 模式”定义偏 微分方程，确定方程类型和方程系数 c,a,f,d，根据具体情况，还可以在不同子区域声明不同系数 . Step 4 “ Mesh 模式”网格化区域 ，可以控制自动生成网格的参数，对生成的网格进行多次细化，使网格分割更细更合理 . Step 5 “ Solve 模式”解偏微分方程，对于椭圆型方程可以激活并控制非线性自适应解题器来处理非线性方程；对于抛物线型方程和双曲型方程，设置初始边界条件后可以求出给定时刻 t 的解；对于特征值问题，可以求出给定区间上的特征值 .求解完成后，可以 返回到 Step 4，对网格进一步细化，进行再次求解
28、 . Step 6 “ View 模式”计算结果的可视化，可以通过设置系统提供的对话框，显示所求的解的表面图、网格图、等高线图和箭头梯形图 .对于抛物线型和双曲线型问题的解还可以进行动画演示 . (2)实例说明用法 求解一个正方形区域上的特征值问题： 12|0u u uu 正方形区域为： 1 1, 1 1.xx (1)使用 PDE 工具箱打开 GUI 求解方程 (2)进入 Draw 模式，绘制一个矩形，然后双击矩形， 在弹出的对话框中设置Left=-1,Bottom=-1,Width=2,Height=2,确认并关闭对话框 (3)进入 Boundary 模式，边界条件采用 Dirichlet
29、条件的默认值 (4)进入 PDE 模式，单击工具栏 PDE 按钮，在弹出的对话框中方程类型选择12 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)单击工具栏的“ =”按钮，开始求解
