求解受限于给出非负解的延迟微分方程 (DDE) 系统

Question

在 MATLAB 中，ode45 有一个名为 NonNegative 的参数，它将解限制为非负数。 They even wrote a paper about how this method works 以及它如何不像将 y_i 变为负数时设置为 0 那样愚蠢，因为这通常不起作用。

现在，MATLAB 也有 dde23 用于求解延迟微分方程，但此积分器没有等效的 NonNegative 参数。

不幸的是，我的任务是向现有 ODE 系统添加延迟，该系统使用 ode45 和启用 NonNegative 解决。

任何想法我应该如何进行？

编辑：

我不确定这是否有帮助，但是...

我系统的 DDE 部分基本上是这样的：

 dx = 1/(1+y*z) - x;
 dy = (y*z)^2/(1+(y*z)^2) - y;
 dz = X - z;

其中X（第三个等式中的大写字母变量）是x 的延迟版本。然后，我通过在x 和z 的方程中添加几个项，然后将组合系统全部集成在一起，将这个 DDE 系统链接到现有的（和更大的）ODE 系统。

Answer 1

您遇到了一个棘手的问题，我不确定是否有一站式解决方案。我很乐意为愿意提供替代答案的任何人提供荣誉。

根据延迟的长度，一种选择是多次运行等式，每次迭代将 x 的旧值传递给最新的更新。

例如，假设您的延迟是一小时。在第一个小时，运行 ode45 并标记 NonNegative。将值与时间参数一起存储到新矩阵中，然后再次运行算法。这次确保添加两个输入参数：旧的解矩阵和旧的时间矩阵

dx = 1/(1+y*z) - x; 
dy = (y*z)^2/(1+(y*z)^2) - y;
tindex = find(told>t,1) -1 % find the upper index which best approximates t
X = xold(tindex) + (xold(tindex+1)-xold(tindex))*(t-told(tindex))/(told(tindex+1)-told(tindex)) % or interpolation method of your choosing
dz = X - z;

现在清洗、冲洗并重复。请注意，X 现在是一个准时间相关项，如来自ode45 的示例 3 所示。

Answer 2

最近我的一些代码遇到了这个问题。 “最简单”的解决方案是执行以下操作，首先，一旦解决方案达到 0，您必须将其保持在 0， 从而改变

dx = 1/(1+y*z) - x; 

to（注意评估 x == 0 案例的位置）

if x > 0 
  dx = 1/(1+y*z) - x; 
else % if x <= 0
  dx = 0;
end

或者可能到（取决于为什么它可能永远不会是 0）

dxTmp = 1/(1+y*z) - x;
if x > 0 
  dx = dxTmp; 
elseif dxTmp > 0
  % x becomes positive again
  dx = dxTmp;
else
  dx = 0;
end

但是请注意，这会在一阶导数中产生跳跃不连续性，当 DDE 求解器到达t - delay 附近的点时，它不会非常有效地求解它，除非它知道该不连续性所在的确切位置（通常你有一个额外的选项告诉 Matlab 它在哪里，但如果你按照下面的步骤，那么就不需要了）。

要确定这种不连续性的位置，您需要使用events 的 DDE 选项（向下滚动到“事件位置属性”，您还可以查看这些examples，其中一个示例实际上处理类似的ODE 中不允许负值的系统 - ODE 和 DDE 的事件几乎相同）。基本上，事件是具有向量输出的 Matlab 函数，向量的每个条目都是对变量的某种或其他评估；在每一步，Matlab 检查其中一个是否等于 0，当其中一个等于 0 时，DDE 停止并返回一个直到该点的解，您必须从该点重新启动 DDE，并将该部分解作为您的历史，即，而不是运行

sol = dde23(ddefun, lags, history, [t0 tEnd], options);

你跑了（注意 sol 和 t0 改变了）

sol = dde23(ddefun, lags, sol, [tCurrent tEnd], options);

在这种情况下，向量的条目之一将是 x（因为您希望 DDE 在 x 等于 0 时停止）。此外，elseif dxTmp <= 0 行代码创建了另一个不连续性，因此当dxTmp 变为0 时，您需要一个事件，即1/(1+y*z) - x 将是向量输出的另一个组成部分。

现在，当您重新启动 ODE 时，Matlab 会自动假定此时存在不连续性，因此您不必担心告诉 Matlab 那里存在不连续性。

下一个问题是，Matlab 从来没有完全正确地实现它，x、y、z 和 X 的值会稍微负一些。因此，如果它会产生问题，您将需要使用

更正 x 的值（以及类似的其他值）

if x < 0
  x = 0;
end

在计算导数之前。但这只会在本地更改 x 的值。因此，您可能还想将 final 解决方案中 x 的所有负值更改为 0。我建议您在将sol 输入到sol = dde23(ddefun, lags, sol, [tCurrent tEnd], options); 之前不要尝试更改它，因为我对此进行了多次尝试，但无法使其正常工作。

Answer 3

有一个非常简单的答案，使用createOptimProblem 来设置优化。您必须使用此方法为每个参数包含边界，但强制参数保持为正值变得微不足道。

详情请看：http://www.mathworks.com/help/gads/createoptimproblem.html