MatlabMatlab教程教程数学科学与技术学院数学科学与技术学院胡金燕胡金燕lionfrtom.com lionfrtom.com应用matlab软件对常微分 方程求解前沿：科学技术和工程中许多问题是用微分方程的形科学技术和工程中许多问题是用微分方程的形 式建立数学模型，因此微分方程的求解有很实式建立数学模型，因此微分方程的求解有很实 际的意义。际的意义。一、常微分方程（组）的符号解一、常微分方程（组）的符号解 二、常微分方程（组）数值解二、常微分方程（组）数值解一、常微分方程（组）的符号解一、常微分方程（组）的符号解函数 dsolve 格式： r = dsolve(eq1,eq2,cond1,cond2,v)说明说明 (1)(1)对给定的常微分方程（组）对给定的常微分方程（组）eq1,eq2,eq1,eq2, 中指定的符号自变量中指定的符号自变量v v，与给定的边界条件，与给定的边界条件 和初始条件和初始条件cond1,cond2,cond1,cond2,. .求符号解（即解求符号解（即解 析解）析解）r r； (2)(2)若没有指定变量若没有指定变量v v，则缺省变量为，则缺省变量为t t；(3)(3)在微分方程（组）的表达式在微分方程（组）的表达式eqeq中，大写中，大写 字母字母D D表示对自变量表示对自变量( (设为设为x)x)的微分算子：的微分算子： D=d/dxD=d/dx，D2=dD2=d2 2/dx/dx2 2，。微分算子。微分算子D D后面的后面的 字母则表示为因变量，即待求解的未知函数字母则表示为因变量，即待求解的未知函数 。初始和边界条件由字符串表示：。初始和边界条件由字符串表示：y(a)=by(a)=b， Dy(c)=dDy(c)=d，D2y(e)=fD2y(e)=f，等等，分别表示，等等，分别表示(4)若边界条件少于方程（组）的阶数 ，则返回的结果r中会出现任意常数 C1，C2，；(5) dsolve命令最多可以接受12个输入参 量（包括方程组与定解条件个数，当然我 们可以做到输入的方程个数多于12个，只 要将多个方程置于一字符串内即可）。(6)若没有给定输出参量，则在命令窗口显 示解列表。若该命令找不到解析解，则返 回一警告信息，同时返回一空的sym对象 。这时，用户可以用命令ode23或 ode45求解方程组的数值解。例1 dsolve(D2y = -a2*y,y(0) = 1,Dy(pi/a) = 0,x) ans= cos(a*x) dsolve(D2y = -a2*y, y(0) = 1,Dy(pi/a) = 0)例2 u,v = dsolve(Du=v,Dv=u)u = C1*exp(-t)+C2*exp(t) V= -C1*exp(-t)+C2*exp(t)二、常微分方程（组）数值解二、常微分方程（组）数值解Matlab专门用于求解常微分方程的函 数，主要采用Runge-Kutta方法： ode23, ode45, ode113, ode15s, ode23s, ode23t, ode23tb二、常微分方程（组）数值解二、常微分方程（组）数值解T,Y = solver(odefun,tspan,y0) T,Y = solver(odefun,tspan,y0,options)T,Y = solver(odefun,tspan,y0,options,p1, p2) 参数说明：参数说明： （1）solver为命令 Ode45,ode23,ode113,ode15s, ode23s,ode23t,ode23tb之一。 （2）odefun 为常微分方程y=f(x,y)， 或为包含一混合矩阵的方程 (x,y)*y=f(x,y). （3）tspan 积分区间（即求解区间）的向 量tspan=t0,tf。要获得问题在其他指定 时间点t0,t1,t2,上的解，则令 tspan=t0,t1,t2,tf （要求是单调的）。参数说明：参数说明： （4）y0 包含初始条件的向量。 （5）options 用命令odeset设置的可选 积分参数. （6）p1,p2, 传递给函数odefun的可选 参数。 在区间在区间tspan=t0,tftspan=t0,tf上，从上，从t0t0到到tftf，用，用 初始条件初始条件y0y0求解显式微分方程求解显式微分方程y y=f(t,y)=f(t,y)。对于标量对于标量t t与列向量与列向量y y，函数，函数f=odefun(t,y)f=odefun(t,y) 必须返回一必须返回一f(t,y)f(t,y)的列向量的列向量f f。解矩阵解矩阵Y Y中的每一行对应于返回的时间列向中的每一行对应于返回的时间列向 量量T T中的一个时间点。中的一个时间点。要获得问题在其他指定时间点要获得问题在其他指定时间点t0,t1,t2,t0,t1,t2, 上的解，则令上的解，则令tspan=t0,t1,t2,tspan=t0,t1,t2,tf,tf（要求是（要求是 单调的）。单调的）。T,Y = solver(odefun,tspan,y0) 用参数用参数optionsoptions（用命令（用命令odesetodeset生成）设置生成）设置 属性（代替了缺省的积分参数），再进行操作属性（代替了缺省的积分参数），再进行操作 。常用的属性包括相对误差值。常用的属性包括相对误差值RelTolRelTol（缺省值（缺省值 为为1e-31e-3）与绝对误差向量）与绝对误差向量AbsTolAbsTol（缺省值为每（缺省值为每 一元素为一元素为1e-61e-6）。）。T,Y = solver(odefun,tspan,y0,options)将参数将参数p1,p2,p3,.p1,p2,p3,.等传递给函数等传递给函数 odefunodefun，再进行计算。若没有参数设置，则，再进行计算。若没有参数设置，则 令令options=options=。T,Y = solver(odefun,tspan,y0,options,p1,p2)求解具体ODE的基本过程：（1）根据问题所属学科中的规律、定律、 公式，用微分方程与初始条件进行描述。 F(y,y,y,y(n),t) = 0 y(0)=y0,y(0)=y1,y(n-1)(0)=yn-1 而y=y;y(1);y(2);,y(m-1)， n与m可以不等求解具体ODE的基本过程：（2 2）运用数学中的变量替换运用数学中的变量替换：y yn n=y=y(n-1)(n-1),y,yn-1n-1=y=y(n-2)(n-2),y,y2 2=y,y=y,y1 1=y=y， 把高阶把高阶（大于大于2 2阶阶）的方程的方程（组组）写成一阶写成一阶 微分方程组微分方程组： （3 3）根据（）根据（1 1）与（）与（2 2）的结果，编写能计算）的结果，编写能计算 导数的导数的M-M-函数文件函数文件odefileodefile。 （4 4）将文件）将文件odefileodefile与初始条件传递给求解与初始条件传递给求解 器器SolverSolver中的一个，运行后就可得到中的一个，运行后就可得到ODEODE 的、在指定时间区间上的解列向量的、在指定时间区间上的解列向量y y（其中包（其中包 含含y y及不同阶的导数）。及不同阶的导数）。不同求解器Solver的特点求解器 SolverODE 类型特点说明ode45非 刚 性一步算法；4，5 阶Runge-Kutta方 程；累计截断误 差达(x)3大部分场合的首 选算法ode23非 刚 性一步算法；2，3 阶Runge-Kutta方 程；累计截断误 差达(x)3使用于精度较低 的情形求解器 SolverOD E类 型特点说明ode113非 刚 性多步法；Adams算 法；高低精度均可 到10-310-6计算时间比ode45 短ode23t适 度 刚 性采用梯形算法适度刚性情形ode15s刚 性多步法；Gears反 向数值微分；精度 中等若ode45失效时， 可尝试使用不同求解器Solver的特点求解器 SolverOD E类 型特点说明ode23s刚 性一步法；2阶 Rosebrock算法； 低精度当精度较低时， 计算时间比 ode15s短ode23tb刚 性梯形算法；低精 度当精度较低时， 计算时间比 ode15s短参数设置属性名取值含义 AbsTol有效值： 正实数或 向量 缺省值： 1e-6绝对误差对应于解向量 中的所有元素；向量则 分别对应于解向量中的 每一分量RelTol有效值： 正实数 缺省值： 1e-3相对误差对应于解向量 中的所有元素。在每步( 第k步)计算过程中，误 差估计为： e(k)1 缺省值：k = 1若k1，则增加每个积分 步中的数据点记录，使 解曲线更加的光滑参数设置属性名取值含义 Jacobian有效值： on、off 缺省值： off若为on时，返回相应的 ode函数的Jacobi矩阵Jpattern有效值： on、off 缺省值： off为on时，返回相应的 ode函数的稀疏Jacobi 矩阵参数设置属性名取值含义 Mass有效值： none、M 、 M(t)、 M(t,y) 缺省值： noneM：不随时间变化的常 数矩阵 M(t)：随时间变化的矩 阵 M(t,y)：随时间、地点变 化的矩阵MaxStep有效值： 正实数 缺省值： tspans/10最大积分步长例3创建函数创建函数function2,function2, 保存在保存在function2.mfunction2.m中中 function f=function2(t,x)function f=function2(t,x) f=-x.2;f=-x.2;在命令窗口中执行在命令窗口中执行 t,x=ode45(function2,0,1,1); t,x=ode45(function2,0,1,1); plot(t,x,-,t,x,o); plot(t,x,-,t,x,o); xlabel(time t0=0,tt=1); xlabel(time t0=0,tt=1); ylabel(x values x(0)=1); ylabel(x values x(0)=1);例4创建函数创建函数function3,function3,保存在保存在function3.mfunction3.m中：中： function f=function3(t,x)function f=function3(t,x) f=x(1)-0.1*x(1)*x(2)+0.01*t;.f=x(1)-0.1*x(1)*x(2)+0.01*t;.-x(2)+0.02*x(1)*x(2)+0.04*t; -x(2)+0.02*x(1)*x(2)+0.04*t;运行命令文件runf3.m t,x=ode45(function3,0,20,30 ;20); plot(t,x); xlabel(time t0=0,tt=20); ylabel(x values x1(0)=30,x2(0)=20);例5创建函数function4，存在function4.m中 function f = function4(t,x) global U f = x(2);U*(1-x(1)2)*x(2)- x(1);再运行命令文件再运行命令文件runf4.mrunf4.m： global Uglobal U U = 7;U = 7; Y0=1;0;