求解器Solver与方程组的关系

森之张卫东

2．求解器Solver与方程组的关系表见表2-3。

表2-3

函数指令		含  义	函  数		含    义
求解器    Solver	ode23	普通2-3阶法解ODE	odefile		包含ODE的文件
	ode23s	低阶法解刚性ODE	选项	odeset	创建、更改Solver选项
	ode23t	解适度刚性ODE		odeget	读取Solver的设置值
	ode23tb	低阶法解刚性ODE	输出	odeplot	ODE的时间序列图
	ode45	普通4-5阶法解ODE		odephas2	ODE的二维相平面图
	ode15s	变阶法解刚性ODE		odephas3	ODE的三维相平面图
	ode113	普通变阶法解ODE		odeprint	在命令窗口输出结果

3．因为没有一种算法可以有效地解决所有的ODE问题，为此，MATLAB提供了多种求解器Solver，对于不同的ODE问题，采用不同的Solver。

表2-4  不同求解器Solver的特点

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

4．在计算过程中，用户可以对求解指令solver中的具体执行参数进行设置（如绝对误差、相对误差、步长等）。

表2-5  Solver中options的属性

属性名	取值	含义
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	有效值：正整数k>1    缺省值：k = 1	若k>1，则增加每个积分步中的数据点记录，使解曲线更加的光滑
Jacobian	有效值：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微分方程file:///C:/Users/lx/AppData/Local/Temp/msohtmlclip1/01/clip_image002.png

y(0)=1，y’(0)=0

令x1=y，x2=dy/dx，则

dx1/dt = x2

dx2/dt = μ(1-x2)-x1

编写函数文件verderpol.m：

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微分方程图