MATLAB 数学应用 微分方程 常微分方程 求解非刚性ODE

本文介绍两个使用 ode45 来求解非刚性常微分方程的示例。MATLAB拥有三个非刚性 ODE 求解器。

• ode45
• ode23
• ode113

对于大多数非刚性问题，ode45 的性能最佳。但对于允许较宽松的误差容限或刚度适中的问题，建议使用 ode23。同样，对于具有严格误差容限的问题，ode113 可能比 ode45 更加高效。

如果非刚性求解器需要很长时间才能解算问题或总是无法完成积分，则该问题可能是刚性问题。有关详细信息，请参阅解算刚性 ODE。

示例：非刚性 van der Pol 方程

van der Pol 方程为二阶 ODE

$$y_1^{\prime \prime }-\mu \left(1-y_1^2\right)y_1^{\prime }+y_1=0,$$

其中 $\mu >0$ 为标量参数。通过执行 $y_1^{\prime }=y_2$ 代换，将此方程重写为一阶 ODE 方程组。生成的一阶 ODE 方程组为

$$\begin{array}{rl}& y^{\prime }=y_2\\ & y_2=\mu (1-y_1^2)y_2-y_1.\end{array}$$

ODE 方程组必须编码为 ODE 求解器能够使用的函数文件。ODE 函数的一般函数签名为

dydt = odefun(t,y)
即，函数必须同时接受 t 和 y 作为输入，即使它没有将 t 用于任何计算时亦如此。

函数文件 vdp1.m 使用 $\mu =1$ 为该 van der Pol 方程编码。变量 $y_1$ 和 $y_2$ 表示为 y(1) 和 y(2)，二元素列向量 dydt 包含 $y_1^{\prime }$ 和 $y_2^{\prime }$ 的表达式。

function dydt = vdp1(t,y)
dydt = [y(2); (1-y(1)^2)*y(2)-y(1)];

使用 ode45 函数、时间区间 [0 20] 和初始值 [2 0] 来解算该 ODE。输出为时间点列向量 t 和解数组 y。y 中的每一行都与 t 的相应行中返回的时间相对应。y 的第一列与 $y_1$ 相对应，第二列与 $y_2$ 相对应。

[t,y] = ode45(@vdp1,[0 20],[2; 0]);

绘制 $y_1$ 和 $y_2$ 的解对 t 的图。

plot(t,y(:,1),'-o',t,y(:,2),'-o')
title('Solution of van der Pol Equation (\mu = 1) using ODE45');
xlabel('Time t');
ylabel('Solution y');
legend('y_1','y_2')

vdpode 函数可求解同一问题，但它接受的是用户指定的 $\mu$ 值。随着 $\mu$ 的增大，van der Pol 方程组将变成刚性。例如，对于值 $\mu =1000$，您需要使用 ode15s 等刚性求解器来求解该方程组。

示例：非刚性欧拉方程

对于专用于非刚性问题的 ODE 求解器，不受外力作用的刚体对应的欧拉方程是其标准测试问题。

这些方程包括

$$\begin{array}{rl}& y_1^{\prime }=y_2{y}_3\\ & y_2^{\prime }=-y_1{y}_3\\ & y_3^{\prime }=-0.51y_1{y}_2.\end{array}$$

函数文件 rigidode 定义此一阶方程组，并在时间区间 [0 12] 上使用初始条件向量 [0; 1; 1]（该向量对应于 $y_1$、$y_2$ 和 $y_3$ 的初始值）对该方程组进行求解。局部函数 f(t,y) 用于编写该方程组的代码。

rigidode 在调用 ode45 时未使用任何输出参数，因此求解器会在每一步之后使用默认的输出函数 odeplot 自动绘制解点。

function rigidode

tspan = [0 12];
y0 = [0; 1; 1];

figure;
ode45(@f,tspan,y0);

% --------------------------------------------------------------------------

function dydt = f(t,y)
dydt = [    y(2)*y(3)
   -y(1)*y(3)
   -0.51*y(1)*y(2) ];

通过调用 rigidode 函数来解算非刚性欧拉方程。

rigidode
title('Solution of Rigid Body w/o External Forces using ODE45')
legend('y_1','y_2','y_3','Location','Best')