24.1 常微分方程 ¶

函数 lsode 可用于求解形式为

dx
-- = f (x, t)
dt

的常微分方程（ODE），使用 Hindmarsh 的 ODE 求解器 LSODE。

[x, istate, msg] = lsode (fcn, x_0, t) ¶

[x, istate, msg] = lsode (fcn, x_0, t, t_crit) ¶

常微分方程（ODE）求解器。

要求解的微分方程组为

dx
-- = f (x, t)
dt

初始条件为

x(t_0) = x_0

解以矩阵 x 的形式返回，每一行对应于向量 t 中的一个元素。t 的第一个元素应为 t_0，且应对应于系统 x_0 的初始状态，因此输出矩阵的第一行就是 x_0。

第一个参数 fcn 是一个字符串、内联函数或函数句柄，用于指定函数 f 的名称，该函数计算方程组的右端向量。该函数必须具有如下形式

xdot = f (x, t)

其中 xdot 和 x 是向量，t 是标量。

如果 fcn 是一个包含两个元素的字符串数组，或者是一个包含两个字符串、内联函数或函数句柄的元胞数组，则第一个元素指定上述函数 f，第二个元素指定一个用于计算 f 的雅可比矩阵的函数。雅可比矩阵函数必须具有如下形式

jac = j (x, t)

其中 jac 是偏导数矩阵

             | df_1  df_1       df_1 |
              | ----  ----  ...  ---- |
              | dx_1  dx_2       dx_N |
              |                       |
              | df_2  df_2       df_2 |
              | ----  ----  ...  ---- |
       df_i   | dx_1  dx_2       dx_N |
jac = ---- = |                       |
       dx_j   |  .    .     .    .    |
              |  .    .      .   .    |
              |  .    .       .  .    |
              |                       |
              | df_M  df_M       df_M |
              | ----  ----  ...  ---- |
              | dx_1  dx_2       dx_N |

第二个参数指定系统的初始状态 x_0。第三个参数是一个向量 t，指定需要求解的时间值。

第四个参数是可选的，可用于指定一组 ODE 求解器不应积分超越的时间点。这对于避免奇点以及导数不连续点处的困难非常有用。

计算成功后，istate 的值将为 2（与 LSODE 的 Fortran 版本一致）。

如果计算不成功，istate 的值将不是 2，而 msg 将包含附加信息。

您可以使用函数 lsode_options 来设置 lsode 的可选参数。

参见 Alan C. Hindmarsh, "ODEPACK, A Systematized Collection of ODE Solvers", Scientific Computing, R. S. Stepleman, editor, 1983. 或 https://computing.llnl.gov/projects/odepack 以获取关于 lsode 内部工作原理的更多信息。

示例：求解范德波尔方程

fvdp = @(y,t) [y(2); (1 - y(1)^2) * y(2) - y(1)];
t = linspace (0, 20, 100);
y = lsode (fvdp, [2; 0], t);

参见： daspk, dassl, dasrt。

lsode_options () ¶

val = lsode_options (opt) ¶

lsode_options (opt, val) ¶

查询或设置 lsode 函数的选项。

当无参数调用时，显示所有可用选项及其当前值。

给定一个参数时，返回选项 opt 的值。

当给定两个参数时，lsode_options 将选项 opt 设置为值 val。

选项包括

"absolute tolerance"

绝对容差。可以是向量或标量。如果是向量，则必须与状态向量的维度匹配。

"relative tolerance"

相对容差参数。与绝对容差不同，此参数只能是标量。

在每个积分步应用的局部误差测试为

  abs (local error in x(i)) <= ...
      rtol * abs (y(i)) + atol(i)

"integration method"

指定用于求解 ODE 系统的积分方法的字符串。有效值为

"adams"

"non-stiff"

不使用雅可比矩阵（即使可用）。

"bdf"

"stiff"

使用刚性向后差分公式（BDF）方法。如果未提供计算雅可比矩阵的函数，lsode 将计算雅可比矩阵的有限差分近似。

"initial step size"

在第一步尝试的步长（默认值自动确定）。

"maximum order"

限制求解方法的最大阶数。如果使用 Adams 方法，此选项必须在 1 到 12 之间。否则，必须在 1 到 5 之间（含）。

"maximum step size"

设置最大步长可以避免跨越非常大的区域（默认值未指定）。

"minimum step size"

允许的最小绝对步长（默认为 0）。

"step limit"

允许的最大步数（默认为 100000）。

"jacobian type"

一个字符串，指定与刚性向后差分公式（BDF）积分方法一起使用的雅可比矩阵类型。有效值为

"full"

默认值。所有偏导数均通过近似或从用户提供的雅可比矩阵函数中获得。

"banded"

仅对对角线以及由选项 "lower jacobian subdiagonals" 和 "upper jacobian subdiagonals" 分别指定的下子对角线和上子对角线数量进行近似或从用户提供的雅可比矩阵函数中使用。用户提供的雅可比矩阵函数可以将所有其他偏导数设置为任意值。

"diagonal"

如果用户提供了雅可比矩阵函数，此设置无效。由 lsode 近似的雅可比矩阵限制在对角线上，其中每个偏导数通过对状态的所有元素同时施加有限变化来计算；如果真实的雅可比矩阵确实始终是对角矩阵，则这与仅对状态的相应元素施加有限变化具有相同的效果，但效率更高。

"lower jacobian subdiagonals"

当选项 "jacobian type" 设置为 "banded" 时使用的下子对角线数量。默认值为 0。

"upper jacobian subdiagonals"

当选项 "jacobian type" 设置为 "banded" 时使用的上子对角线数量。默认值为 0。

以下是一个使用 lsode 求解三个微分方程组的示例。给定函数

## oregonator 微分方程
function xdot = f (x, t)

  xdot = zeros (3,1);

  xdot(1) = 77.27 * (x(2) - x(1)*x(2) + x(1) ...
            - 8.375e-06*x(1)^2);
  xdot(2) = (x(3) - x(1)*x(2) - x(2)) / 77.27;
  xdot(3) = 0.161*(x(1) - x(3));

endfunction

以及初始条件 x0 = [ 4; 1.1; 4 ]，可以使用以下命令对方程组进行积分

t = linspace (0, 500, 1000);

y = lsode ("f", x0, t);

如果您尝试此操作，您会看到结果值在 t = 0 到 5 之间以及 t ≈ 305 附近发生剧烈变化。一组更高效的输出点可能是

t = [0, logspace(-1, log10(303), 150), ...
        logspace(log10(304), log10(500), 150)];

上述微分方程对应的 m 文件包含在 Octave 发行版的示例目录中，名称为 oregonator.m。