当前位置：文江博客话题详情

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

发布于 2024-11-28 12:29:41 字数 1202 浏览 3 评论 0原文

在 MATLAB 中，ode45 有一个名为非负这将解决方案限制为非负。他们甚至写了一篇论文，介绍这种方法的工作原理以及它为何不是那么愚蠢的事情每当 y_i 变为负数时，只需将 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 系统。系统全部在一起。

原文

In MATLAB, ode45 has a parameter called NonNegative which constrains the solutions to be nonnegative. They even wrote a paper about how this method works and how it's not something as stupid as just setting y_i to 0 whenever it becomes negative, as that won't generally work.

Now, MATLAB also has dde23 for solving delay differential equations, but there is no equivalent NonNegative parameter for this integrator.

Unfortunately, I am tasked with adding a delay to an existing ODE system which is solved using ode45 with NonNegative enabled.

Any ideas how I should proceed?

EDIT:

I'm not sure if this is helpful, but...

The DDE part of my system basically looks like:

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

where X (the capital letter variable in the third equation) is the delayed version of x. Then, I'm linking this DDE system to an existing (and larger) ODE system by adding a couple terms to the equations for x and z, and then integrating the combined system all together.

分享到QQ

分享到微博

如果你对这篇内容有疑问，欢迎到本站社区发帖提问参与讨论，获取更多帮助，或者扫码二维码加入 Web 技术交流群。

发布评论

需要登录才能够评论，你可以免费注册一个本站的账号。

宛菡 2024-12-05 12:29:41

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

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

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

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 所示。

You have a tough problem and I'm not sure if there is a one-step solution. I'd be more then glad to provide kudos to anyone willing to provide an alternative answer.

Depending on the length of the delay, one option would be to run the equation several times, with each iteration passing the old values of x to latest update.

For instance, say your delay is one hour. In the first hour, run ode45 with the NonNegative flagged. Store the Value into a New matrix along with the Time parameter, and run the algorithm again. This time make sure you add two input parameters: your old solution matrix and the old time matrix

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;

Now wash, rinse, and repeat. Note that X is now a quasi-time-dependant term as seen in example 3 from ode45.

回复收藏 0 原文

梦初启 2024-12-05 12:29:41

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

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

为（注意 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 - 延迟的点，它不能非常有效地解决它，除非它知道这个不连续性的确切位置（通常你使用一个额外的选项来告诉 Matlab它在哪里，但是如果如果您按照以下步骤操作，则不需要）。

要确定此不连续性的位置，您需要使用事件（向下滚动到“活动地点属性”，您还可以查看这些示例，其中一个示例实际上涉及类似的系统，其中 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 的值会略有负值。因此，如果它会产生问题，您需要

if x < 0
  x = 0;
end

在计算导数之前更正 x 的值（以及类似的其他值）。但这只会在本地更改 x 的值。因此，您可能还想将最终解决方案中x的所有负值更改为0。我建议您在将 sol 输入到 sol = dde23(ddefun, lags, sol, [tCurrent tEnd], options); 之前不要尝试更改 sol，因为我做了一个尝试了很多次，但我无法让它发挥作用。

I recently had this problem with some code of mine. The 'simplest' solution is to do the following, firstly once the solution reaches 0, you have to keep it at 0,
thus change

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

to (notice where the x == 0 case is evaluated)

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

or maybe to (depending on why it may never be 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

Note however that this creates a jump discontinuity in the first derivative, which when the DDE solver reaches the point of having t - delay near this point, it does not solve it very efficiently unless it knows the exact spot where this discontinuity is (usually you have use an extra option of telling Matlab where it is, but if you follow the steps below, then that will not be needed).

To determine the place of this discontinuity you need to use the DDE option of events (scroll down to 'Event Location Property', you can also look at these examples, one of those examples actually deals with a similar system where negative values are not allowed in an ODE - events for ODEs and DDEs are almost identical). Basically an event is a Matlab function with a vector output, each one of the entries of the vector is some or other evaluation of your variables; at each step Matlab checks if one of them eqauls 0, when one of them equal 0 the DDE stops and returns a solution up to that point, from which you must restart the DDE with that partial solution as your history, i.e. instead of running

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

you run (notice sol and t0 changed)

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

In this case one of the vector's entries is going to be x (as you want the DDE to stop when x equals 0). Also the line of code elseif dxTmp <= 0 creates another discontinuity, thus you need an event for when dxTmp becomes 0, i.e. 1/(1+y*z) - x will be another component of the vector output.

Now when you restart the ODE, Matlab automatically assumes that there is a discontinuity at that point, hence you do not have to worry about telling Matlab that there is one there.

The next problem is that Matlab never quite achieves it right, the values of x, y, z and X will be slightly negative. Thus if it is going to create a problem, you will want to correct the value of x (and the other values similarly) with

if x < 0
  x = 0;
end

before calculating the derivatives. But this only changes the value of x locally. So you might want to change all negative values of x in the final solution to 0 as well. I suggest that you do not try to change sol before inputting it into sol = dde23(ddefun, lags, sol, [tCurrent tEnd], options); as I made a number of attempts at that and I could not get it to work.

回复收藏 0 原文