求解非刚性 ODE

打开脚本

本页包含两个使用 ode45 来求解非刚性常微分方程的示例。MATLAB® 提供几个非刚性 ODE 求解器。

  • ode45

  • ode23

  • ode78

  • ode89

  • ode113

对于大多数非刚性问题，ode45 的性能最佳。但对于允许较宽松的误差容限或刚度适中的问题，建议使用 ode23。同样，对于具有更严格误差容限的问题，或当计算 ODE 函数的计算成本很高时，ode113 可能比 ode45 更高效。ode78 和 ode89 是高阶求解器，在精度对稳定性至关重要的长时积分中表现出色。

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

示例：非刚性范德波尔方程

范德波尔方程为二阶 ODE

$$y''_1 - \mu \left( 1 - y_1^2\right) y'_1+y_1=0,$$

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

$$
\begin{array}{cl}
y'_1 &= 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$ 为该范德波尔方程编码。变量 $y_1$ 和 $y_2$ 表示为 y(1) 和 y(2)，二元素列向量 dydt 包含 $y'_1$ 和 $y'_2$ 的表达式。

function dydt = vdp1(t,y)
%VDP1  Evaluate the van der Pol ODEs for mu = 1
%
%   See also ODE113, ODE23, ODE45.

%   Jacek Kierzenka and Lawrence F. Shampine
%   Copyright 1984-2014 The MathWorks, Inc.

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$ 的增大，范德波尔方程组将变成刚性。例如，对于值 $\mu = 1000$，您需要使用 ode15s 等刚性求解器来求解该方程组。

示例：非刚性欧拉方程

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

这些方程包括

$$ \begin{array}{cl} y'_1 &= y_2y_3 \\ y'_2 &= -y_1y_3 \\ y'_3 &=
-0.51y_1y_2. \end{array}$$

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

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

function rigidode
%RIGIDODE  Euler equations of a rigid body without external forces.
%   A standard test problem for non-stiff solvers proposed by Krogh.  The
%   analytical solutions are Jacobian elliptic functions, accessible in
%   MATLAB.  The interval here is about 1.5 periods; it is that for which
%   solutions are plotted on p. 243 of Shampine and Gordon.
%
%   L. F. Shampine and M. K. Gordon, Computer Solution of Ordinary
%   Differential Equations, W.H. Freeman & Co., 1975.
%
%   See also ODE45, ODE23, ODE113, FUNCTION_HANDLE.

%   Mark W. Reichelt and Lawrence F. Shampine, 3-23-94, 4-19-94
%   Copyright 1984-2014 The MathWorks, Inc.

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

% solve the problem using ODE45
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')

另请参阅

| | | |

主题

外部网站