ddesd giving NaN for all state variables beyond time step 1.

The following is my code for a state-dependent delay equation I am working on.

% Parameter set and intialization
%(Q=20, a=5, k=1, ell=10)
tstar = 60.2762; %tstar for b=.03
tspan = [tstar 5000];

options = ddeset(RelTol=1e-12); %Without this extreme tolerance, ddesd never runs at all
sol = ddesd(@ddefun, @dely, @history, tspan, options); %call ddesd

% Print graph
axes
hold on
plot(sol.x,sol.y);
legend('C','n1','m1','n2','m2','tau')
hold off
xlim([0 400])
xlabel('Time (s)')
ylabel('Total Fission Rate (nM/s)')

%ODE which defines behavior before critical time tstar
function dydt = odefun(t,y)
b = .03;
k = 1;

dydt = [-k*y(1)^2 - k*y(1)*y(2) + b*y(3);
    k*y(1)^2 - b*y(2);
    k*y(1)*(y(1) + y(2)) - b*y(3);
    ];
end

function dydt = ddefun(t,y,Z) % DDE being solved
b = .03;
k = 1;
a = 5;
ell = 10;
dydt = [a*y(5) - k*y(1) - k*y(1)*(y(2) + y(4)) + b*(y(3) + y(4));
    -k*y(1)*(Z(1,1)*exp(-b*y(6)) - y(1)) - b*y(2);
    -k*y(1)*(ell*Z(1,1)*exp(-b*y(6)) - y(1) - y(2)) - b*y(3);
    k*y(1)*Z(1,1)*exp(-b*y(6)) - a*y(4);
    k*y(1)*(ell*Z(1,1)*exp(-b*y(6)) + y(4)) - a*y(5);
    1 - y(1)/Z(1,1)
    ];

end
%-------------------------------------------
function d = dely(t,y) % delay for y
d = t - y(6);
end
%-------------------------------------------
function v = history(t) % history function for t < t0
tstar = 60.2762; %tstar for b=.03

preTstarSol = ode45(@odefun, [0 tstar], [20 0 0]); %Here the ODE is called to create the history of C, n1, and m1

IC1 = @(t) interp1(preTstarSol.x, preTstarSol.y(1,:), t, 'linear');
IC2 = @(t) interp1(preTstarSol.x, preTstarSol.y(2,:), t, 'linear');
IC3 = @(t) interp1(preTstarSol.x, preTstarSol.y(3,:), t, 'linear');

v = [IC1(t);
    IC2(t);
    IC3(t);
    0;
    0;
    tstar]; %The delay is initialized as "tstar" meaning in the first time step the delay looks back to time 0

end
%-------------------------------------------

As it is written, the solver finds all the state variables to be NaN after time step 1, and prints an empty graph. You may note the rather extreme RelTol=1e-12. However, without a change to the tolerance, nothing happens at all. That is, if you delete "options" from inside ddesd, the code will never finish running and even when you click "STOP" there is no "sol" initialized at all, which implies the code never even runs the first time step of ddesd.

What I can't understand is that a tiny tweak to parameter b and the corresponding tweak to tstar: tstar1 = 37.8275; %tstar for b = .05, it runs just fine and produces the damped oscillations I expect from the system.

Any ideas what could be causing the issue?

Answer 1

[Torsten]

The history function has to cover values for t that are negative (because y6 becomes large and you have to supply y1(t-y6)). But you only supply history data for y1(t), y2(t) and y3(t) between 0 and tstar for t. This means that IC1(t), IC2(t) and IC3(t) will return NaN. Something like

t = 0:0.1:1;
y = t.^2;
interp1(t,y,-2,'linear')

ans = NaN

[Sam Chak]

Hi @Torsten

I am checking my understanding. The OP intends to solve the DDE system of six states from s to 5000 s. The state‑dependent delay in this system is . Six solution histories are provided via the history() function so ddesd solver can access values for times preceding the initial integration point s. The first three histories are time‑dependent (particularly for ) and the final three histories are constant.

At the 1st integration step , ddesd can evaluate the delay term at the initial time by referencing at , i.e. . However, at the 2nd integration step , since , the state increases and it is possible that , so refers to at a negative time , which is undefined and resulting in getting NaN values.

tstar = 60.2762;
tspan = linspace(0, tstar, 603);
v     = history(tspan)

v = 3×603

20.0000 3.9256 1.1259 0.3836 0.1662 0.0895 0.0675 0.0597 0.0570 0.0562 0.0560 0.0561 0.0562 0.0563 0.0565 0.0567 0.0568 0.0570 0.0572 0.0573 0.0575 0.0577 0.0579 0.0580 0.0582 0.0584 0.0585 0.0587 0.0589 0.0590 0 10.2768 10.8226 10.8462 10.8205 10.7897 10.7579 10.7260 10.6942 10.6624 10.6308 10.5992 10.5677 10.5364 10.5051 10.4740 10.4429 10.4118 10.3810 10.3502 10.3195 10.2888 10.2583 10.2279 10.1974 10.1673 10.1371 10.1070 10.0771 10.0473 0 16.0744 18.8741 19.6164 19.8338 19.9105 19.9325 19.9403 19.9430 19.9438 19.9440 19.9439 19.9438 19.9437 19.9435 19.9433 19.9432 19.9430 19.9428 19.9427 19.9425 19.9423 19.9421 19.9420 19.9418 19.9416 19.9415 19.9413 19.9411 19.9410

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

% initial value of y1 in the DDE system at t0 = tstar
v(1,end)

ans = 0.2160

function dydt = odefun(t,y)
    b   = .03;
    k   = 1;
    dydt= [-k*y(1)^2 - k*y(1)*y(2) + b*y(3);
            k*y(1)^2 - b*y(2);
            k*y(1)*(y(1) + y(2)) - b*y(3)];
end

% solution function for t < t0
function v = history(t)
    tstar = 60.2762;
    
    preTstarSol = ode45(@odefun, [0 tstar], [20; 0; 0]);
    
    IC1 = @(t) interp1(preTstarSol.x, preTstarSol.y(1,:), t, 'linear');
    IC2 = @(t) interp1(preTstarSol.x, preTstarSol.y(2,:), t, 'linear');
    IC3 = @(t) interp1(preTstarSol.x, preTstarSol.y(3,:), t, 'linear');
    
    v = [IC1(t);
        IC2(t);
        IC3(t)];
end

[Kitrick]

Hello both. Thank you for engaging with my question. If I am understanding what you both are saying, I should try to expand my history function into negative time. Here is my attempt, and unfortunately I still get NaN for all state variables for all time steps beyond 1.

%(Q=20, a=5, k=1, ell=10)
tstar = 60.2762; %tstar for b=.03
tspan = [tstar 5000];

preTstarSol = ode45(@odefun, [0 tstar], [20 0 0]);

options = ddeset(RelTol=1e-12);
sol = ddesd(@ddefun, @dely, @history, tspan, options);

axes
hold on
plot(sol.x,sol.y);
legend('C','n1','m1','n2','m2','tau')
hold off
xlim([0 400])
xlabel('Time (s)')
ylabel('Total Fission Rate (nM/s)')

plot(preTstarSol.x,preTstarSol.y);
legend('C','n1','m1');

function dydt = odefun(t,y)
b = .03;
k = 1;

dydt = [-k*y(1)^2 - k*y(1)*y(2) + b*y(3);
    k*y(1)^2 - b*y(2);
    k*y(1)*(y(1) + y(2)) - b*y(3);
    ];

end

function dydt = ddefun(t,y,Z) % equation being solved
b = .03;
k = 1;
a = 5;
ell = 10;
dydt = [a*y(5) - k*y(1) - k*y(1)*(y(2) + y(4)) + b*(y(3) + y(4));
    -k*y(1)*(Z(1,1)*exp(-b*y(6)) - y(1)) - b*y(2);
    -k*y(1)*(ell*Z(1,1)*exp(-b*y(6)) - y(1) - y(2)) - b*y(3);
    k*y(1)*Z(1,1)*exp(-b*y(6)) - a*y(4);
    k*y(1)*(ell*Z(1,1)*exp(-b*y(6)) + y(4)) - a*y(5);
    1 - y(1)/Z(1,1)
    ];

end
%-------------------------------------------
function d = dely(t,y) % delay for y
d = t - y(6);
end
%-------------------------------------------
function v = history(t) % history function for t < t0
tstar = 60.2762; %tstar for b=.03

preTstarSol = ode45(@odefun, [0 tstar], [20 0 0]);

% Here is the change I made, having each initial condition function be
% constant at it's time=0 value for all time between -10 and 0.
IC1 = @(t) interp1([-10, preTstarSol.x], [preTstarSol.y(1,1), preTstarSol.y(1,:)], t, 'linear');
IC2 = @(t) interp1([-10, preTstarSol.x], [preTstarSol.y(2,1), preTstarSol.y(2,:)], t, 'linear');
IC3 = @(t) interp1([-10, preTstarSol.x], [preTstarSol.y(3,1), preTstarSol.y(3,:)], t, 'linear');

v = [IC1(t);
    IC2(t);
    IC3(t);
    0;
    0;
    tstar];

end
%-------------------------------------------

Is there an issue with my implementation of history into negative time? Or does this show something else is wrong?

[Torsten]

Resetting the "ddeset" option of "RelTol=1e-12" makes negative time values in "history" disappear. I nonetheless used an if-statement to avoid possible NaN values returned from "interp1" for the case that parameters might be changed.

Thus your original code would have worked (at least up to t = 800) if you hadn't requested such a stringent precision.

%(Q=20, a=5, k=1, ell=10)
tstar = 60.2762; %tstar for b=.03
tspan = [tstar 800];

preTstarSol = ode45(@odefun, [0 tstar], [20 0 0]);

options = [];%ddeset(RelTol=1e-12);
sol = ddesd(@ddefun, @dely, @history, tspan, options);

axes
hold on
plot(sol.x,sol.y);
legend('C','n1','m1','n2','m2','tau')
hold off
%xlim([0 400])
xlabel('Time (s)')
ylabel('Total Fission Rate (nM/s)')

plot(preTstarSol.x,preTstarSol.y);
legend('C','n1','m1');

function dydt = odefun(t,y)
b = .03;
k = 1;

dydt = [-k*y(1)^2 - k*y(1)*y(2) + b*y(3);
    k*y(1)^2 - b*y(2);
    k*y(1)*(y(1) + y(2)) - b*y(3);
    ];

end

function dydt = ddefun(t,y,Z) % equation being solved
b = .03;
k = 1;
a = 5;
ell = 10;
dydt = [a*y(5) - k*y(1) - k*y(1)*(y(2) + y(4)) + b*(y(3) + y(4));
    -k*y(1)*(Z(1,1)*exp(-b*y(6)) - y(1)) - b*y(2);
    -k*y(1)*(ell*Z(1,1)*exp(-b*y(6)) - y(1) - y(2)) - b*y(3);
    k*y(1)*Z(1,1)*exp(-b*y(6)) - a*y(4);
    k*y(1)*(ell*Z(1,1)*exp(-b*y(6)) + y(4)) - a*y(5);
    1 - y(1)/Z(1,1)
    ];

end
%-------------------------------------------
function d = dely(t,y) % delay for y
d = t - y(6);
end
%-------------------------------------------
function v = history(t) % history function for t < t0
tstar = 60.2762; %tstar for b=.03
if t < 0
  v = [20;
       0;
       0;
       0;
       0;
       tstar]; 
  return
end

preTstarSol = ode45(@odefun, [0 tstar], [20 0 0]);

% Here is the change I made, having each initial condition function be
% constant at it's time=0 value for all time between -10 and 0.
IC1 = interp1(preTstarSol.x, preTstarSol.y(1,:), t, 'linear');
IC2 = interp1(preTstarSol.x, preTstarSol.y(2,:), t, 'linear');
IC3 = interp1(preTstarSol.x, preTstarSol.y(3,:), t, 'linear');

v = [IC1;
     IC2;
     IC3;
     0;
     0;
     tstar];
end
%-------------------------------------------

[Kitrick]

Thank you again for your help. This has fixed the issue I was having. It is with some embarassment that I must admit the real issue for why the solution was not behaving the way I was expecting is that the line:

dydt = [a*y(5) - k*y(1) - k*y(1)*(y(2) + y(4)) + b*(y(3) + y(4));

should actually be

dydt = [a*y(5) - k*y(1)^2 - k*y(1)*(y(2) + y(4)) + b*y(3);

Another great reminder to double, triple, quadruple check that you copied the system from your notes to the computer correctly.