6

在 MATLAB 中，ode45有一个名为的参数NonNegative将解约束为非负数。他们甚至写了一篇关于这种方法如何工作的论文，以及它如何不像将 y_i 设置为 0 那样愚蠢，因为它通常不会起作用。

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

不幸的是，我的任务是向现有的 ODE 系统添加延迟，该系统可以使用ode45with NonNegativeenabled 解决。

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

编辑：

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

我系统的 DDE 部分基本上如下所示：

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

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

4

3 回答 3

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 它在哪里的额外选项，但如果您按照以下步骤操作，则不需要）。

要确定这种不连续性的位置，您需要使用事件的 DDE 选项（向下滚动到“事件位置属性”，您还可以查看这些示例，其中一个示例实际上处理的是不允许负值的类似系统在 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);

在这种情况下，向量的条目之一将是（因为您希望 DDE 在等于 0x时停止）。x代码行还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本地的值。因此，您可能还希望将最终x解决方案中的所有负值更改为 0。我建议您在输入之前不要尝试更改，因为我对此进行了多次尝试，但无法使其正常工作。solsol = dde23(ddefun, lags, sol, [tCurrent tEnd], options);

于 2013-12-02T02:12:09.203 回答
2

你有一个棘手的问题，我不确定是否有一步解决方案。我会更乐意为任何愿意提供替代答案的人提供荣誉。

根据延迟的长度，一种选择是多次运行等式，每次迭代将 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 所示。

于 2011-08-08T21:07:29.213 回答
-2

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

详细信息：http: //www.mathworks.com/help/gads/createoptimproblem.html

于 2015-04-09T16:26:57.557 回答