matlab解常微分方程的特别方法

Matlab教程数学科学与技术学院胡金燕[e-mail]应用matlab软件对常微分方程求解前沿：科学技术和工程中许多问题是用微分方程的形式建立数学模型，因此微分方程的求解有很实际的意义。一、常微分方程（组）的符号解二、常微分方程（组）数值解一、常微分方程（组）的符号解函数dsolve格式：r=dsolve('eq1,eq2,…’,'cond1,cond2,…','v')说明(1)对给定的常微分方程（组）eq1,eq2,…中指定的符号自变量v，与给定的边界条件和初始条件cond1,cond2,….求符号解（即解析解）r；(2)若没有指定变量v，则缺省变量为t；(3)在微分方程（组）的表达式eq中，大写字母D表示对自变量(设为x)的微分算子：D=d/dx，D2=d2/dx2，…。微分算子D后面的字母则表示为因变量，即待求解的未知函数。初始和边界条件由字符串表示：y(a)=b，Dy(c)=d，D2y(e)=f，等等，分别表示b)x(yaxd)x(ycxf)x(yex(4)若边界条件少于方程（组）的阶数，则返回的结果r中会出现任意常数C1，C2，…；(5)dsolve命令最多可以接受12个输入参量（包括方程组与定解条件个数，当然我们可以做到输入的方程个数多于12个，只要将多个方程置于一字符串内即可）。(6)若没有给定输出参量，则在命令窗口显示解列表。若该命令找不到解析解，则返回一警告信息，同时返回一空的sym对象。这时，用户可以用命令ode23或ode45求解方程组的数值解。例1dsolve('D2y=-a^2*y,y(0)=1,Dy(pi/a)=0','x')0)/(1)0(2ayyyayans=cos(a*x)dsolve(‘D2y=-a^2*y’,‘y(0)=1,Dy(pi/a)=0’)例2[u,v]=dsolve('Du=v,Dv=u')uvvuu=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,tf]上，从t0到tf，用初始条件y0求解显式微分方程y’=f(t,y)。对于标量t与列向量y，函数f=odefun(t,y)必须返回一f(t,y)的列向量f。解矩阵Y中的每一行对应于返回的时间列向量T中的一个时间点。要获得问题在其他指定时间点t0,t1,t2,…上的解，则令tspan=[t0,t1,t2,…,tf]（要求是单调的）。[T,Y]=solver(odefun,tspan,y0)用参数options（用命令odeset生成）设置属性（代替了缺省的积分参数），再进行操作。常用的属性包括相对误差值RelTol（缺省值为1e-3）与绝对误差向量AbsTol（缺省值为每一元素为1e-6）。[T,Y]=solver(odefun,tspan,y0,options)将参数p1,p2,p3,..等传递给函数odefun，再进行计算。若没有参数设置，则令options=[]。[T,Y]=solver(odefun,tspan,y0,options,p1,p2…)求解具体ODE的基本过程：（1）根据问题所属学科中的规律、定律、公式，用微分方程与初始条件进行描述。F(y,y’,y’’,…,y(n),t)=0y(0)=y0,y’(0)=y1,…,y(n-1)(0)=yn-1而y=[y;y(1);y(2);…,y(m-1)]，n与m可以不等求解具体ODE的基本过程：（2）运用数学中的变量替换：yn=y(n-1),yn-1=y(n-2),…,y2=y‘,y1=y，把高阶（大于2阶）的方程（组）写成一阶微分方程组：),(),(),(2121ytfytfytfyyyynn（3）根据（1）与（2）的结果，编写能计算导数的M-函数文件odefile。（4）将文件odefile与初始条件传递给求解器Solver中的一个，运行后就可得到ODE的、在指定时间区间上的解列向量y（其中包含y及不同阶的导数）。nnyyyyyyy10210)0()0()0(不同求解器Solver的特点求解器SolverODE类型特点说明ode45非刚性一步算法；4，5阶Runge-Kutta方程；累计截断误差达(△x)3大部分场合的首选算法ode23非刚性一步算法；2，3阶Runge-Kutta方程；累计截断误差达(△x)3使用于精度较低的情形求解器SolverODE类型特点说明ode113非刚性多步法；Adams算法；高低精度均可到10-3～10-6计算时间比ode45短ode23t适度刚性采用梯形算法适度刚性情形ode15s刚性多步法；Gear’s反向数值微分；精度中等若ode45失效时，可尝试使用不同求解器Solver的特点求解器SolverODE类型特点说明ode23s刚性一步法；2阶Rosebrock算法；低精度当精度较低时，计算时间比ode15s短ode23tb刚性梯形算法；低精度当精度较低时，计算时间比ode15s短参数设置属性名取值含义AbsTol有效值：正实数或向量缺省值：1e-6绝对误差对应于解向量中的所有元素；向量则分别对应于解向量中的每一分量RelTol有效值：正实数缺省值：1e-3相对误差对应于解向量中的所有元素。在每步(第k步)计算过程中，误差估计为：e(k)=max(RelTol*abs(y(k)),AbsTol(k))参数设置属性名取值含义NormControl有效值：on、off缺省值：off为‘on’时，控制解向量范数的相对误差，使每步计算中，满足：norm(e)=max(RelTol*norm(y),AbsTol)Events有效值：on、off为‘on’时，返回相应的事件记录参数设置属性名取值含义OutputFcn有效值：odeplot、odephas2、odephas3、odeprint缺省值：odeplot若无输出参量，则solver将执行下面操作之一：画出解向量中各元素随时间的变化；画出解向量中前两个分量构成的相平面图；画出解向量中前三个分量构成的三维相空间图；随计算过程，显示解向量参数设置属性名取值含义OutputSel有效值：正整数向量缺省值：[]若不使用缺省设置，则OutputFcn所表现的是那些正整数指定的解向量中的分量的曲线或数据。若为缺省值时，则缺省地按上面情形进行操作Refine有效值：正整数k1缺省值：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1)0('2xxx创建函数function2,保存在function2.m中functionf=function2(t,x)f=-x.^2;在命令窗口中执行[t,x]=ode45('function2',[0,1],1);plot(t,x,'-',t,x,'o');xlabel('timet0=0,tt=1');ylabel('xvaluesx(0)=1');00.20.40.60.810.50.550.60.650.70.750.80.850.90.951timet0=0,tt=1xvaluesx(0)=1例420)0(30)0(04.002.0'01.01.0'2121222111xxtxxxxtxxxx创建函数function3,保存在function3.m中：functionf=function3(t,x)f=[x(1)-0.1*x(1)*x(2)+0.01*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('timet0=0,tt=20');ylabel('xvaluesx1(0)=30,x2(0)=20');05101520020406080100120timet0=0,tt=20xvaluesx1(0)=30,x2(0)=20例50)0(',1)0(01)1(222yydtdyydtyd1)1(,2212222121xxdtyddtdxxdtdydtdxdtdyxyx，则令创建函数function4，存在function4.m中functionf=function4(t,x)globalUf=[x(2);U*(1-x(1)^2)*x(2)-x(1)];再运行命令文件runf4.m：globalUU=7;Y0=[1;0];[t,x]=ode45('function4',0,40,Y0);x1=x(:,1);x2=x(:,2);plot(t,x1,t,x2)0510152025303540-15-10-5051015