有没有办法在 MATLAB 中输出 ode 求解器的中间计算？

Question

我有一个刚性 ODE 系统，我用 ode15 求解。时间导数的类型为：

yt = A*y - u.*y;

其中

yt

、

y

和

u

是列向量，

A

是稀疏矩阵。每次评估时都必须进行

u

的计算，并且计算量很大（u 是非线性的，它是刚度的原因）。我不想预先计算

u

，因为我想保留求解器稍后选择任何所需时间的自由。然而，我感兴趣的不是解 y(t)，而是 u(t).*y(t)，它已经被计算来计算

yt

。如果我得到一个结构体

sol

的解决方案，然后我可以用

deval

进行评估，我可以轻松计算 y(t)，但必须重新计算 u(t)。我想知道是否有一种方法可以在使用求解器时存储

u.*y

的结果，这样我就可以像调用

deval(y,t)

一样轻松地调用它。我发现的最接近的替代方案是依靠 deval 也提供时间导数，所以然后做：

[y,yt] = deval(sol,t); %t is a row vector
uy     = A*y - yt;     %uy = u.*y

这并不是非常低效，但在已经计算过的操作中浪费了一些计算。

知道如何做到这一点吗？

Answer 1

据我研究了几个月前遇到的完全相同的问题，根据内置函数这是不可能的。有一些解决方法可以保存中间计算（例如使用持久变量或全局变量），但它们会产生更多问题。

MATLAB 论坛上有很多讨论。 这个答案可能特别有用。其他评论可以在这里或这里找到。

使用 MATLAB 内置

ode

函数保存中间计算的问题是，它们是 自适应步长，并且每步要调用多次要积分的函数。 如果将

t

作为具有 2 个以上元素的向量传递，MATLAB 将按提供的时间步长返回解。然而，这并不意味着它会在这些时间步长上计算解决方案 - 它只是从其内置解决方案中进行插值。

最好的选择就是你提到的，在解决颂歌后计算你想要的结果，这意味着重新计算

u

。如果你正确地向量化你的函数，这一步应该不会有时间问题。

Answer 2

您可以“破解”ODE 求解器并将控制变量

u

包含在状态向量中。然后你的 ODE 系统就变成了一个（简单的、显式的）DAE 系统

y' = f(t,y,u)
0  = u - g(t,x)

现在唯一的“技巧”是，您还需要在选项中传递一个质量矩阵，该矩阵是一个对角矩阵，其条目

1

用于微分方程，

0

用于代数方程。检查文档您使用的求解器是否支持 DAE，所有隐式方法都应提供此功能。

可能会有一些效率损失，因为新方程也会影响步长控制，但我预计不会有很大差异。

在结果中，控件的值将存储在状态向量的最后一个分量中，因此提取它们只是对该向量进行简单的切片。

Answer 3

我已经按照 Lutz Lehmann 在此线程中建议的破解 ODE 求解器的方法取得了成功。但请注意，此方法需要奇异质量矩阵，这意味着您必须使用刚性求解器，例如 ode15s，按照最初的问题。

非刚性求解器无法获得奇异质量矩阵，例如 ode45，因此此方法不适用于非刚性求解器。