《matlab解常微分方程的特别方法课件》由会员分享，可在线阅读，更多相关《matlab解常微分方程的特别方法课件（35页珍藏版）》请在金锄头文库上搜索。
1、MatlabMatlab教程教程教程教程数学科学与技术学院数学科学与技术学院数学科学与技术学院数学科学与技术学院 胡金燕胡金燕胡金燕胡金燕 应用matlab软件对常微分方程求解前沿：前沿：科学技术和工程中许多问题是用微分方程的形科学技术和工程中许多问题是用微分方程的形式建立数学模型，因此微分方程的求解有很实式建立数学模型，因此微分方程的求解有很实际的意义。际的意义。一、常微分方程（组）的符号解一、常微分方程（组）的符号解二、常微分方程（组）数值解二、常微分方程（组）数值解一、常微分方程（组）的符号解一、常微分方程（组）的符号解函数函数 dsolve 格式：格式：r = dsolve(eq1,e
2、q2,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
3、=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个输入参个输入参量（包括方程组与定解条件个数，当然我量
4、（包括方程组与定解条件个数，当然我们可以做到输入的方程个数多于们可以做到输入的方程个数多于12个，只个，只要将多个方程置于一字符串内即可）。要将多个方程置于一字符串内即可）。(6)若没有给定输出参量，则在命令窗口显若没有给定输出参量，则在命令窗口显示解列表。若该命令找不到解析解，则返示解列表。若该命令找不到解析解，则返回一警告信息，同时返回一空的回一警告信息，同时返回一空的sym对象。对象。这时，用户可以用命令这时，用户可以用命令ode23或或ode45求解方程组的数值解。求解方程组的数值解。例例1 dsolve(D2y = -a2*y,y(0) = 1,Dy(pi/a) = 0,x)ans=
5、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 = sol
6、ver(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。要获
7、得问题在其他指定。要获得问题在其他指定时间点时间点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)。 对于标量对于标
8、量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) 用参数用参数optionsoptio
9、ns（用命令（用命令odesetodeset生成）设置生成）设置属性（代替了缺省的积分参数），再进行操作。属性（代替了缺省的积分参数），再进行操作。常用的属性包括相对误差值常用的属性包括相对误差值RelTolRelTol（缺省值为（缺省值为1e-31e-3）与绝对误差向量）与绝对误差向量AbsTolAbsTol（缺省值为每一（缺省值为每一元素为元素为1e-61e-6）。）。T,Y = solver(odefun,tspan,y0,options) 将参数将参数p1,p2,p3,.等传递给函数等传递给函数odefun，再进行计算。若没有参数设置，则，再进行计算。若没有参数设置，则令令option
10、s=。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 2）运用数学中的变量替换：）运用数学中的变量替换：yn=y(n-1),yn-1=y(n-2),y2=y
11、,y1=y，把高阶（大于把高阶（大于2 2阶）的方程（组）写成一阶阶）的方程（组）写成一阶微分方程组：微分方程组： （3）根据（）根据（1）与（）与（2）的结果，编写能计算）的结果，编写能计算导数的导数的M-函数文件函数文件odefile。（4）将文件）将文件odefile与初始条件传递给求解与初始条件传递给求解器器Solver中的一个，运行后就可得到中的一个，运行后就可得到ODE的、在指定时间区间上的解列向量的、在指定时间区间上的解列向量y（其中包（其中包含含y及不同阶的导数）。及不同阶的导数）。不同求解器不同求解器Solver的特点的特点求解器求解器SolverODE类型类型特点特点说明说
12、明ode45ode45非非刚刚性性一步算法；一步算法；4 4，5 5阶阶Runge-KuttaRunge-Kutta方方程；累计截断误程；累计截断误差达差达(x)(x)3 3大部分场合的首大部分场合的首选算法选算法ode23ode23非非刚刚性性一步算法；一步算法；2 2，3 3阶阶Runge-KuttaRunge-Kutta方方程；累计截断误程；累计截断误差达差达(x)(x)3 3使用于精度较低使用于精度较低的情形的情形求解器求解器SolverODE类类型型特点特点说明说明ode113ode113非非刚刚性性多步法；多步法；AdamsAdams算算法；高低精度均可法；高低精度均可到到1010
13、-3-31010-6-6计算时间比计算时间比ode45ode45短短ode23tode23t适适度度刚刚性性采用梯形算法采用梯形算法适度刚性情形适度刚性情形ode15sode15s刚刚性性多步法；多步法；GearGears s反向数值微分；精反向数值微分；精度中等度中等若若ode45ode45失效时，失效时，可尝试使用可尝试使用不同求解器不同求解器Solver的特点的特点求解器求解器SolverODE类类型型特点特点说明说明ode23sode23s刚刚性性一步法；一步法；2 2阶阶RosebrockRosebrock算法；算法；低精度低精度当精度较低时，当精度较低时，计算时间比计算时间比ode
14、15sode15s短短ode23tbode23tb刚刚性性梯形算法；低精梯形算法；低精度度当精度较低时，当精度较低时，计算时间比计算时间比ode15sode15s短短参数设置参数设置属性名属性名取值取值含义含义AbsTol有效值：有效值：正实数或正实数或向量向量缺省值：缺省值：1e-6绝对误差对应于解向量绝对误差对应于解向量中的所有元素；向量则中的所有元素；向量则分别对应于解向量中的分别对应于解向量中的每一分量每一分量RelTol有效值：有效值：正实数正实数缺省值：缺省值：1e-3相对误差对应于解向量相对误差对应于解向量中的所有元素。在每步中的所有元素。在每步(第第k步步)计算过程中，误计算过
15、程中，误差估计为：差估计为：e(k)=max(RelTol*abs(y(k),AbsTol(k)参数设置参数设置属性名属性名取值取值含义含义NormControl有效值：有效值：on、off缺省值：缺省值：off为为on时，控制解向时，控制解向量范数的相对误差，使量范数的相对误差，使每步计算中，满足：每步计算中，满足：norm(e)1缺省值：缺省值：k = 1若若k1，则增加每个积分，则增加每个积分步中的数据点记录，使步中的数据点记录，使解曲线更加的光滑解曲线更加的光滑参数设置参数设置属性名属性名取值取值含义含义Jacobian有效值：有效值：on、off缺省值：缺省值：off若为若为on时，
16、返回相时，返回相应的应的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创建
17、函数创建函数function2, 保存在保存在function2.m中中function f=function2(t,x)f=-x.2;在命令窗口中执行在命令窗口中执行 t,x=ode45(function2,0,1,1); plot(t,x,-,t,x,o); xlabel(time t0=0,tt=1); ylabel(x values x(0)=1);例例4创建函数创建函数function3,保存在保存在function3.m中：中：function f=function3(t,x)f=x(1)-0.1*x(1)*x(2)+0.01*t;. -x(2)+0.02*x(1)*x(2)+0.
18、04*t;运行命令文件运行命令文件runf3.mt,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 Uf = x(2);U*(1-x(1)2)*x(2)-x(1);再运行命令文件再运行命令文件runf4.m：global UU = 7;Y0=1;0;t,x = ode45(function4,0,40,Y0);x1=x(:,1);x2=x(:,2);plot(t,x1,t,x2)
