常微分方程数值解2.4

《常微分方程数值解2.4》由会员分享，可在线阅读，更多相关《常微分方程数值解2.4（3页珍藏版）》请在金锄头文库上搜索。

1、 1 2.4 常微分方程数值解常微分方程数值解 函数 ode45、ode23、ode113、ode15s、ode23s、ode23t、ode23tb 功能 常微分方程（ODE）组初值问题的数值解 参数说明： solver 为命令 ode45、ode23,ode113,ode15s,ode23s,ode23t,ode23tb 之一。 Odefun 为显式常微分方程 y=f(t,y)，或为包含一混合矩阵的方程 M(t,y)*y=f(t,y)。命令 ode23 只能求解常数混合矩阵的问题；命令 ode23t 与 ode15s 可以求解奇异矩 阵的问题。 Tspan 积分区间（即求解区间）的向量 ts

2、pan=t0,tf。要获得问题在其他指定时间点 t0,t1,t2,上的解，则令 tspan=t0,t1,t2,tf（要求是单调的） 。 Y0 包含初始条件的向量。 Options 用命令 odeset 设置的可选积分参数。 P1,p2, 传递给函数 odefun 的可选参数。 格式 T,Y = solver(odefun,tspan,y0) %在区间 tspan=t0,tf上，从 t0 到 tf，用初始条 件y0求解显式微分方程y=f(t,y)。 对于标量t与列向量y， 函数f=odefun(t,y) 必须返回一 f(t,y)的列向量 f。解矩阵 Y 中的每一行对应于返回的时间列 向量 T 中

3、的一个时间点。 要获得问题在其他指定时间点 t0,t1,t2,上的解， 则令 tspan=t0,t1,t2,tf（要求是单调的） 。 T,Y = solver(odefun,tspan,y0,options) %用参数 options（用命令 odeset 生成） 设置的属性（代替了缺省的积分参数） ，再进行操作。常用的属性包括相 对误差值 RelTol（缺省值为 1e-3）与绝对误差向量 AbsTol（缺省值为每 一元素为 1e-6） 。 T,Y =solver(odefun,tspan,y0,options,p1,p2) 将参数 p1,p2,p3,等传递给函数 odefun，再进行计算。若

4、没有参数设置，则令 options=。 1求解具体 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 可以不等 （2）运用数学中的变量替换：yn=y(n-1),yn-1=y(n-2),y2=y1=y，把高阶（大于 2 阶）的方 程（组）写成一阶微分方程组： = = )y, t (f )y, t (f )y, t (f y y y y n 2 1 n 2 1 M M ， = = n 1 0

5、 n 2 1 0 y y y )0(y )0(y )0(y y MM （3）根据（1）与（2）的结果，编写能计算导数的 M-函数文件 odefile。 （4）将文件 odefile 与初始条件传递给求解器 Solver 中的一个，运行后就可得到 ODE 的、在指定时间区间上的解列向量 y（其中包含 y 及不同阶的导数） 。 2求解器 Solver 与方程组的关系表见表 2-3。 表2-3 函数指令 含 义 函 数 含 义 ode23 普通2-3阶法解ODE odefile 包含ODE的文件 ode23s 低阶法解刚性ODE odeset 创建、更改Solver选项 ode23t 解适度刚性OD

6、E 选项 odeget 读取Solver的设置值 ode23tb 低阶法解刚性ODE odeplot ODE的时间序列图 求解器 Solver ode45 普通4-5阶法解ODE 输出 odephas2 ODE的二维相平面图 2 ode15s 变阶法解刚性ODE odephas3 ODE的三维相平面图 ode113 普通变阶法解ODE odeprint 在命令窗口输出结果 3因为没有一种算法可以有效地解决所有的 ODE 问题，为此，MATLAB 提供了多种 求解器 Solver，对于不同的 ODE 问题，采用不同的 Solver。 表2-4 不同求解器Solver的特点 求解器Solver O

7、DE类型 特点 说明 ode45 非刚性 一步算法；4，5阶Runge-Kutta方 程；累计截断误差达(x)3 大部分场合的首选算法 ode23 非刚性 一步算法；2，3阶Runge-Kutta方 程；累计截断误差达(x)3 使用于精度较低的情形 ode113 非刚性 多步法；Adams算法；高低精度均 可到10-310-6 计算时间比ode45短 ode23t 适度刚性 采用梯形算法 适度刚性情形 ode15s 刚性 多步法；Gears反向数值微分；精 度中等 若ode45失效时，可尝试使 用 ode23s 刚性 一步法；2阶Rosebrock算法；低 精度 当精度较低时，计算时间比 od

8、e15s短 ode23tb 刚性 梯形算法；低精度 当精度较低时，计算时间比 ode15s短 4在计算过程中，用户可以对求解指令 solver 中的具体执行参数进行设置（如绝对误 差、相对误差、步长等） 。 表2-5 Solver中options的属性 属性名 取值 含义 AbsTol 有效值：正实数或向量 缺省值：1e-6 绝对误差对应于解向量中的所有元素；向量则分别对应于 解向量中的每一分量 RelTol 有效值：正实数 缺省值：1e-3 相对误差对应于解向量中的所有元素。 在每步(第k步)计算 过程中，误差估计为： e(k)1，则增加每个积分步中的数据点记录，使解曲线更 加的光滑 Jac

9、obian 有效值：on、off 缺省值：off 若为on时，返回相应的ode函数的Jacobi矩阵 Jpattern 有效值：on、off 缺省值：off 为on时，返回相应的ode函数的稀疏Jacobi矩阵 Mass 有效值：none、M、 M(t)、M(t,y) 缺省值：none M：不随时间变化的常数矩阵 M(t)：随时间变化的矩阵 M(t,y)：随时间、地点变化的矩阵 MaxStep 有效值：正实数 缺省值：tspans/10 最大积分步长 例 2-45 求解描述振荡器的经典的 Ver der Pol 微分方程01 dt dy )y1 ( dt yd 2 2 2 =+ y(0)=1，y(0)=0 令x1=y，x2=dy/dx，则 dx1/dt = x2 dx2/dt = (1-x2)-x1 编写函数文件 verderpol.m： 3 function xprime = verderpol(t,x) global MU xprime = x(2);MU*(1-x(1)2)*x(2)-x(1); 再在命令窗口中执行： global MU MU = 7; Y0=1;0 t,x = ode45(verderpol,0,40,Y0); x1=x(:,1);x2=x(:,2); plot(t,x1,t,x2) 图形结果为图 2-20。 图2-20 Ver der Pol微分方程图