Solving a delay differential equation (DDE) system constrained to give nonnegative solutions

Question

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.

Answer 1

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.

Answer 2

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.

Answer 3

There's a much simpler answer, use createOptimProblem to set up an optimization. You have to include bounds for every parameter using this method, but it becomes trivial to force a parameter to remain positive.

Details here: http://www.mathworks.com/help/gads/createoptimproblem.html