24.2 微分代数方程 ¶

函数 daspk 可用于求解如下形式的 DAE（微分代数方程）：

0 = f (x-dot, x, t),    x(t=0) = x_0, x-dot(t=0) = x-dot_0

其中 x-dot 是 x 的导数。该方程使用 Petzold 的 DAE 求解器 DASPK 求解。

[x, xdot, istate, msg] = daspk (fcn, x_0, xdot_0, t, t_crit) ¶

求解一组微分代数方程。

daspk 求解方程组

0 = f (x, xdot, t)

满足

x(t_0) = x_0, xdot(t_0) = xdot_0

解以矩阵 x 和 xdot 的形式返回，结果矩阵中的每一行对应向量 t 中的一个元素。t 的第一个元素应为 t_0，对应于系统的初始状态 x_0 及其导数 xdot_0，因此输出 x 的第一行是 x_0，输出 xdot 的第一行是 xdot_0。

第一个参数 fcn 是一个字符串、内联函数或函数句柄，用于命名调用以计算方程组残差向量的函数 f。它必须具有如下形式：

res = f (x, xdot, t)

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

如果 fcn 是一个二元素字符串数组或由字符串、内联函数或函数句柄组成的二元素元胞数组，则第一个元素命名上述函数 f，第二个元素命名一个用于计算修正雅可比矩阵的函数：

      df       df
jac = -- + c ------
      dx     d xdot

修正雅可比函数必须具有如下形式：

jac = j (x, xdot, t, c)

daspk 的第二和第三个参数指定状态及其导数的初始条件，第四个参数指定所需解的输出时间向量，包括与初始条件对应的时间。

初始状态和导数的集合并不严格要求必须一致。如果它们不一致，则必须使用 daspk_options 函数提供附加信息，以便 daspk 能够计算一致的起点。

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

成功计算后，istate 的值将大于零（与 DASPK 的 Fortran 版本一致）。

如果计算不成功，istate 的值将小于零，并且 msg 将包含附加信息。

您可以使用函数 daspk_options 为 daspk 设置可选参数。

另请参阅: dassl.

daspk_options () ¶

val = daspk_options (opt) ¶

daspk_options (opt, val) ¶

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

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

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

当用两个参数调用时，daspk_options 将选项 opt 设置为值 val。

选项包括：

"absolute tolerance"

绝对容差。可以是向量或标量。如果是向量，它必须与状态向量的维度匹配，并且相对容差也必须是相同长度的向量。

"relative tolerance"

相对容差。可以是向量或标量。如果是向量，它必须与状态向量的维度匹配，并且绝对容差也必须是相同长度的向量。

每个积分步骤中应用的局部误差测试为：

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

"compute consistent initial condition"

用 ‘Y_d’ 表示状态向量中的微分变量，用 ‘Y_a’ 表示代数变量，ddaspk 可以解决以下两个初始化问题之一：

1. 给定 Y_d，计算 Y_a 和 Y’_d
2. 给定 Y’，计算 Y。

在两种情况下，都需要输入给定分量的初值，并且必须同时提供未知分量的初始猜测值。将此选项设置为 1 可解决第一个问题，设置为 2 可解决第二个问题（默认值为 0，因此您必须提供一组一致的初始条件）。

如果此选项设置为非零值，则还必须设置 "algebraic variables" 选项，以声明问题中的哪些变量是代数变量。

"use initial condition heuristics"

设置为非零值以使用下面描述的初始条件启发式选项。

"initial condition heuristics"

一个由以下参数组成的向量，可用于控制初始条件计算。

MXNIT

牛顿迭代的最大次数（默认值为 5）。

MXNJ

雅可比矩阵求值的最大次数（默认值为 6）。

MXNH

如果 "compute consistent initial condition" 选项已设置为 1，则尝试的人工步长参数值的最大数量（默认值为 5）。

注意，如果 "compute consistent initial condition" 选项已设置为 1，则允许的牛顿迭代最大总数为 MXNIT*MXNJ*MXNH；如果设置为 2，则为 MXNIT*MXNJ。

LSOFF

设置为非零值以禁用线搜索算法（默认值为 0）。

STPTOL

线搜索算法中的最小缩放步长（默认值为 eps^(2/3)）。

EPINIT

牛顿迭代收敛性测试中的摆动因子。该测试应用于残差向量（预乘以近似雅可比矩阵）。为了收敛，该向量的加权 RMS 范数（按误差权重缩放）必须小于 EPINIT*EPCON，其中 EPCON = 0.33 是在时间步中使用的类似测试常数。默认值为 EPINIT = 0.01。

"print initial condition info"

将此选项设置为非零值以显示初始条件计算的详细信息（默认值为 0）。

"exclude algebraic variables from error test"

设置为非零值以从误差测试中排除代数变量。您还必须设置 "algebraic variables" 选项，以声明问题中的哪些变量是代数变量（默认值为 0）。

"algebraic variables"

一个与状态向量长度相同的向量。非零元素表示状态向量的对应元素是代数变量（即其导数在方程组中未显式出现）。

"compute consistent initial condition" 和 "exclude algebraic variables from error test" 选项需要此选项。

"enforce inequality constraints"

设置为以下值之一，以强制实施由 "inequality constraint types" 选项指定的不等式约束（默认值为 0）。

1. 仅在初始条件计算中进行约束检查。
2. 在积分过程中强制进行约束检查。
3. 同时强制实施选项 1 和 2。

"inequality constraint types"

一个与状态向量长度相同的向量，用于指定不等式约束的类型。向量的每个元素对应状态的一个元素，应分配以下代码之一：

-2

小于零。

-1

小于或等于零。

0

无约束。

1

大于或等于零。

2

大于零。

此选项仅在 "enforce inequality constraints" 选项非零时有效。

"initial step size"

微分代数问题有时在第一步可能会出现严重的缩放困难。如果您对问题的缩放有较多了解，可以通过指定初始步长来帮助缓解此问题（默认值为自动计算）。

"maximum order"

限制求解方法的最大阶数。此选项必须介于 1 和 5 之间（包括边界，默认值为 5）。

"maximum step size"

设置最大步长将避免越过非常大的区域（默认未指定）。

Octave 还包含 DASSL（DASPK 的早期版本）和 DASRT，后者可用于求解带约束（停止条件）的 DAE。

[x, xdot, istate, msg] = dassl (fcn, x_0, xdot_0, t, t_crit) ¶

求解一组微分代数方程。

dassl 求解方程组

0 = f (x, xdot, t)

满足

x(t_0) = x_0, xdot(t_0) = xdot_0

解以矩阵 x 和 xdot 的形式返回，结果矩阵中的每一行对应向量 t 中的一个元素。t 的第一个元素应为 t_0，对应于系统的初始状态 x_0 及其导数 xdot_0，因此输出 x 的第一行是 x_0，输出 xdot 的第一行是 xdot_0。

第一个参数 fcn 是一个字符串、内联函数或函数句柄，用于命名调用以计算方程组残差向量的函数 f。它必须具有如下形式：

res = f (x, xdot, t)

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

如果 fcn 是一个二元素字符串数组或由字符串、内联函数或函数句柄组成的二元素元胞数组，则第一个元素命名上述函数 f，第二个元素命名一个用于计算修正雅可比矩阵的函数：

      df       df
jac = -- + c ------
      dx     d xdot

修正雅可比函数必须具有如下形式：

jac = j (x, xdot, t, c)

dassl 的第二和第三个参数指定状态及其导数的初始条件，第四个参数指定所需解的输出时间向量，包括与初始条件对应的时间。

初始状态和导数的集合并不严格要求必须一致。然而在实践中，DASSL 不太擅长为您确定一致集，因此最好确保初始值使函数求值为零。

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

成功计算后，istate 的值将大于零（与 DASSL 的 Fortran 版本一致）。

如果计算不成功，istate 的值将小于零，并且 msg 将包含附加信息。

您可以使用函数 dassl_options 为 dassl 设置可选参数。

另请参阅: daspk, dasrt, lsode.

dassl_options () ¶

val = dassl_options (opt) ¶

dassl_options (opt, val) ¶

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

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

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

当用两个参数调用时，dassl_options 将选项 opt 设置为值 val。

选项包括：

"absolute tolerance"

绝对容差。可以是向量或标量。如果是向量，它必须与状态向量的维度匹配，并且相对容差也必须是相同长度的向量。

"relative tolerance"

相对容差。可以是向量或标量。如果是向量，它必须与状态向量的维度匹配，并且绝对容差也必须是相同长度的向量。

每个积分步骤中应用的局部误差测试为：

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

"compute consistent initial condition"

如果非零，dassl 将尝试计算一组一致的初始条件。这通常不太可靠，因此最好提供一致的集合并将此选项保持为零。

"enforce nonnegativity constraints"

如果您知道方程的解始终为非负，设置此参数为非零值可能会有所帮助。但最好先尝试将此选项保持为零，如果效果不佳再将其设置为非零值。

"initial step size"

微分代数问题有时在第一步可能会出现严重的缩放困难。如果您对问题的缩放有较多了解，可以通过指定初始步长来帮助缓解此问题。

"maximum order"

限制求解方法的最大阶数。此选项必须介于 1 和 5 之间（包括边界）。

"maximum step size"

设置最大步长将避免越过非常大的区域（默认未指定）。

"step limit"

单次调用底层 Fortran 代码时尝试的最大积分步数。

[x, xdot, t_out, istat, msg] = dasrt (fcn, g, x_0, xdot_0, t) ¶

… = dasrt (fcn, g, x_0, xdot_0, t, t_crit) ¶

… = dasrt (fcn, x_0, xdot_0, t) ¶

… = dasrt (fcn, x_0, xdot_0, t, t_crit) ¶

求解一组微分代数方程。

dasrt 求解方程组

0 = f (x, xdot, t)

满足

x(t_0) = x_0, xdot(t_0) = xdot_0

具有函数停止准则（求根功能）。

解以矩阵 x 和 xdot 的形式返回，结果矩阵中的每一行对应向量 t_out 中的一个元素。t 的第一个元素应为 t_0，对应于系统的初始状态 x_0 及其导数 xdot_0，因此输出 x 的第一行是 x_0，输出 xdot 的第一行是 xdot_0。

向量 t 提供了积分长度的上限。如果满足停止条件，向量 t_out 将比 t 短，并且 t_out 的最后一个元素将是满足停止条件的点，可能不与向量 t 中的任何元素对应。

第一个参数 fcn 是一个字符串、内联函数或函数句柄，用于命名调用以计算方程组残差向量的函数 f。它必须具有如下形式：

res = f (x, xdot, t)

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

如果 fcn 是一个二元素字符串数组或由字符串、内联函数或函数句柄组成的二元素元胞数组，则第一个元素命名上述函数 f，第二个元素命名一个用于计算修正雅可比矩阵的函数：

      df       df
jac = -- + c ------
      dx     d xdot

修正雅可比函数必须具有如下形式：

jac = j (x, xdot, t, c)

可选的第二个参数命名一个函数，该函数定义在积分过程中需要求根的约束函数。此函数必须具有如下形式：

g_out = g (x, t)

并返回约束函数值的向量。如果任何约束函数的值改变符号，DASRT 将尝试在符号改变点停止积分。

如果省略了约束函数的名称，dasrt 解决与 daspk 或 dassl 相同的问题。

注意，由于舍入误差和积分误差导致的约束函数数值误差，DASRT 可能返回虚假根，或在两个或多个几乎相等的 T 值处返回相同的根。如果怀疑存在此类虚假根，用户应考虑在约束函数的求值中使用更小的误差容限或更高的精度。

如果某个约束函数的根定义了问题的结束，DASRT 的输入仍应允许积分稍微超过该根点，以便 DASRT 可以通过插值定位根。

dasrt 的第三和第四个参数指定状态及其导数的初始条件，第四个参数指定所需解的输出时间向量，包括与初始条件对应的时间。

初始状态和导数的集合并不严格要求必须一致。然而在实践中，DASSL 不太擅长为您确定一致集，因此最好确保初始值使函数求值为零。

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

成功计算后，istate 的值将大于零（与 DASSL 的 Fortran 版本一致）。

如果计算不成功，istate 的值将小于零，并且 msg 将包含附加信息。

您可以使用函数 dasrt_options 为 dasrt 设置可选参数。

另请参阅: dasrt_options, daspk, dasrt, lsode.

dasrt_options () ¶

val = dasrt_options (opt) ¶

dasrt_options (opt, val) ¶

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

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

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

当用两个参数调用时，dasrt_options 将选项 opt 设置为值 val。

选项包括：

"absolute tolerance"

绝对容差。可以是向量或标量。如果是向量，它必须与状态向量的维度匹配，并且相对容差也必须是相同长度的向量。

"relative tolerance"

相对容差。可以是向量或标量。如果是向量，它必须与状态向量的维度匹配，并且绝对容差也必须是相同长度的向量。

每个积分步骤中应用的局部误差测试为：

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

"initial step size"

微分代数问题有时在第一步可能会出现严重的缩放困难。如果您对问题的缩放有较多了解，可以通过指定初始步长来帮助缓解此问题。

"maximum order"

限制求解方法的最大阶数。此选项必须介于 1 和 5 之间（包括边界）。

"maximum step size"

设置最大步长将避免越过非常大的区域。

"step limit"

单次调用底层 Fortran 代码时尝试的最大积分步数。

请参阅 K. E. Brenan 等人所著的 Numerical Solution of Initial-Value Problems in Differential-Algebraic Equations（北荷兰出版社， 1989 年，DOI: https://doi.org/10.1137/1.9781611971224）， 以获取有关 DASSL 实现的更多信息。