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

Question

Kitrick 2025-11-6，0:25

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2181072-ddesd-giving-nan-for-all-state-variables-beyond-time-step-1

编辑： Kitrick about 15 hours 前

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?

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Torsten 2025-11-6，2:07

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2181072-ddesd-giving-nan-for-all-state-variables-beyond-time-step-1#answer_1571695

编辑：Torsten 2025-11-6，10:09

在 MATLAB Online 中打开

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

4 个评论
显示 2更早的评论隐藏 2更早的评论

Sam Chak 12 minutes 前

在 MATLAB Online 中打开

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 about 20 hours 前

编辑：Kitrick about 20 hours 前

在 MATLAB Online 中打开

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 about 19 hours 前

编辑：Torsten about 18 hours 前

在 MATLAB Online 中打开

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;

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;

tstar];

end

%-------------------------------------------

Kitrick about 17 hours 前

编辑：Kitrick about 15 hours 前

在 MATLAB Online 中打开

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.

请先登录，再进行评论。

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

4 个评论
显示 2更早的评论隐藏 2更早的评论

更多回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

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

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

4 个评论 显示 2更早的评论隐藏 2更早的评论

更多回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

4 个评论
显示 2更早的评论隐藏 2更早的评论