ans =
Help implementing 4th-order Runge–Kutta for ODEs with Dirac impulses (discontinuous system)
Show older comments
Hello
13 Comments
In what you posted I cannot find an obvious mistake. Maybe a unit error ? It's always problematic not to use SI units for all quantities involved (e.g. mg/kg for ui - this implies that the unit for Mi is kg and for U is mg) .
It would help if you posted your complete code together with a link to the article you are referring to.
Walter Roberson
on 1 Dec 2025
That code outline looks reasonable to me.
Paul
on 2 Dec 2025
IMO, trying to anything like this inside the ODE function is not what should be done.
Is mass_fun1(t) continuous at each t_i ?
Torsten
on 2 Dec 2025
Your approach splitting the interval of integration at the times where the dosing of U takes place is the correct way to numerically handle the discontinuities for U at specified times.
A numerical integrator will never exactly capture the time instant where
any(kroneckerDelta(sym(t - tDose(j)), sym(0)))
becomes 1. Thus your U won't have jumps at the times where dosing takes place - you will simply solve
dUdt = -kYU * U
If you post code and a link to the article, I will take a look if I can find discrepancies.
ode45() uses variable time steps. It does not increment time by any fixed value such as 0.01; instead it takes larger and larger steps until prediction fails, in which case it remains where it is and takes increasingly short steps until prediction succeeds again... and then goes back to taking larger and larger steps. Accordingly, it is not very likely to happen to match any given time exactly -- and without exact matches, kroneckerDelta is not going to work.
Even if, hypothetically, it incremented time by fixed steps of (say) 0.001 then there would be the question of round-off error in accumulating all of those 0.001 . 0.001 has no exact representation in binary floating point; the closest to it is
fprintf('%.999g\n', 0.001)
add several of those together and the fraction is going to accumulate... leading to the distinct possibility that by (for example) time 2.5 seconds that the accumulated time does not match 2.5 exactly ... and again, kroneckerDelta would only work for bit-for-bit matches.
We have thus established that kroneckerDelta is unlikely to work consistently no matter whether the timesteps are fixed or variable.
Hi @Júlia
Is it acceptable to model a Dirac-type impulse as a Gaussian whose width gets very small while keeping its area equal to one?
% parameters
M1.kYU = 1; % dU/dt is self-stabilizing, where U converges to zero, if kYU > 0
M1.kEY = 1; % Since U → 0 as t → ∞, if kEY > 0, then dY/dt is stable and Y → 0
M1.kZY = 1; % The sign does not affect the stability, after all Y → 0 as t → ∞
M1.kEZ = 1; % Since Y → 0 as t → ∞, if kEZ > 0, then dZ/dt is stable and Z → 0
M1.kX = 1;
M1.kEX = 2; % Since Y → 0 as t → ∞, if kEX >kX, then dX/dt is stable and X → 0
M1.etaEY = 1; % The sign does not affect the stability, after all Y → 0 as t → ∞
M1.lambda = 1; % will affect how fast the term exp(- λ·Z) converges to 0 as t → ∞
% Mouse dynamics
function dydt = fM1_fun(t,y,M1)
% approximation of a Dirac-type impulse as a Gaussian function
amp = 7978.85; % amplitude of the Gaussian
sig = 1e-4; % width of the Gaussian
ctr = 0; % center of the Gaussian
dose= amp*gaussmf(t, [sig, ctr]); % the area of this function is one
% ODEs
dU = -M1.kYU * y(1) + dose;
dY = M1.kYU * y(1) - M1.kEY * y(2);
dZ = M1.kZY * y(2) - M1.kEZ * y(3);
dX = M1.kX * y(4) - M1.kEX * y(4) - M1.etaEY * y(2) .* y(4) .* exp(-M1.lambda * y(3));
dydt= [dU; dY; dZ; dX];
end
% solve the system
U0 = 1.0; % because the area of the impulse-like signal is 1, U0 will appear at 2
Y0 = 1/3;
Z0 = - 1/3;
X0 = - 1.0;
[t, y] = ode45(@(t, y) fM1_fun(t,y,M1), [0 10], [U0; Y0; Z0; X0]);
plot(t, y), grid on
legend('U', 'Y', 'Z', 'X')
xlabel('t')
ylabel('Amplitude')
@Júlia, if you're comparing against a published paper then including a citation to the journal that contains the paper (issue, title, etc.) in a comment would allow anyone with access to that journal to see the difference between your plots and the paper's plots.
At first glance, I'm wondering about the values in your M1 struct. You wrote "even though the ODEs, parameters, dose times, and mass values were reproduced as closely as possible from the article." -- how precisely were those parameters given in the article? These values are pretty close and both round to 5.5 to one decimal place:
lambdaRange = [5.4501; 5.5499]
diff(lambdaRange)
round(lambdaRange, 1)
But if you were to use them as the input to the exp function, the values you get are relatively far apart.
format longg
E = exp(lambdaRange)
diff(E)
Sometimes small differences in inputs can lead to larger differences in outputs.
If it's convenient for you, could you share your email so I can send you the document and the full code?
Sorry, but I don't want to share my email address.
But it should be no problem to solve for U, Y and Z analytically. The equation for X is nonlinear - thus more complicated.
%syms U(t) Y(t) Z(t) X(t)
%syms kYU kEY kZY kEZ kX kEX eta lambda real
syms U(t) Y(t) Z(t)
syms kYU kEY kZY kEZ real
eqn1 = diff(U,t) == -kYU*U;
eqn2 = diff(Y,t) == kYU*U - kEY*Y;
eqn3 = diff(Z,t) == kZY*Y - kEZ*Z;
%eqn4 = diff(X,t) == kX*X - kEX*X - eta*Y*X*exp(-lambda*Z)
sol = dsolve([eqn1,eqn2,eqn3]);
%sol = dsolve([eqn1,eqn2,eqn3,eqn4])
sol.U
sol.Y
sol.Z
%sol.X
@Torsten, the changes in the solutions are likely undetectable to the naked eye when the impulsive dose is absent and the initial value for U is adjusted to 2. However, when the impulsive-like dose is administered to Mouse U at 10 seconds after the states have settled, it appears that the states remain unaffected. @Júlia, are these 'restful' behaviors expected (see Case 1 below)?
Since the solutions for the first three linear differential equations can be determined, can the solution for 
be expressed in this manner?
be expressed in this manner?
By the way, @Steven Lord made an excellent observation regarding the sensitivity of the exponential function. As a reviewer, I sometimes attempt to replicate the results reported in the manuscript. If there are significant deviations, and when I request clarification from the authors, they often acknowledge that some parameter values in the Table have been rounded or truncated due to space constraints in the two-column journal format. The authors should include a footnote, marked with an asterisk (*), to clarify the true values used.
In Case 2, the impulse function is approximated using a relatively larger sigma value in the Gaussian function while ensuring that the area under the curve remains equal to one. The expected 'jerk' results are obtained. However, the state X does not appear to be affected because the time-varying term 
does not contain state U where the impulsive-like dose is administered, and the state Y already converges to a significantly small value.
does not contain state U where the impulsive-like dose is administered, and the state Y already converges to a significantly small value.% parameters
M1.kYU = 1; % dU/dt is self-stabilizing, where U converges to zero, if kYU > 0
M1.kEY = 1; % Since U → 0 as t → ∞, if kEY > 0, then dY/dt is stable and Y → 0
M1.kZY = 1; % The sign does not affect the stability, after all Y → 0 as t → ∞
M1.kEZ = 1; % Since Y → 0 as t → ∞, if kEZ > 0, then dZ/dt is stable and Z → 0
M1.kX = 1;
M1.kEX = 2; % Since Y → 0 as t → ∞, if kEX >kX, then dX/dt is stable and X → 0
M1.etaEY = 1; % The sign does not affect the stability, after all Y → 0 as t → ∞
M1.lambda = 1; % will affect how fast the term exp(- λ·Z) converges to 0 as t → ∞
% Mice dynamics
function [dydt, dose] = fM1_fun(t, y, M1, i) % with the i-Switch
% approximation of a Dirac-type impulse δ(t - 10) as a Gaussian function
if i == 1
sig = 1e-3; % relatively low value of sigma
else
sig = 4e-3; % sigma: width of the Gaussian (don't go lower than 4e-3)!
end
% amp = 2/(sig*sqrt(2*pi)); % amplitude of the Gaussian (if δ(t) is injected at t = 0)
amp = 1/(sig*sqrt(2*pi)); % amplitude of the Gaussian (if δ(t) is time-shifted)
ctr = 10; % center of the Gaussian
dose= amp*gaussmf(t, [sig, ctr]); % the area under this function is one
% ODEs
dU = -M1.kYU * y(1) + dose;
dY = M1.kYU * y(1) - M1.kEY * y(2);
dZ = M1.kZY * y(2) - M1.kEZ * y(3);
dX = M1.kX * y(4) - M1.kEX * y(4) - M1.etaEY * y(2) .* y(4) .* exp(-M1.lambda * y(3));
dydt= [dU; dY; dZ; dX];
end
% solve the Mice system
U0 = 1.0;
Y0 = 1/3;
Z0 = - 1/3;
X0 = - 1.0;
tspan = 0:1e-5:20;
options = odeset(RelTol=1e-8, AbsTol=1e-10);
% Case 1
i = 1;
[t, y] = ode45(@(t, y) fM1_fun(t, y, M1, i), tspan, [U0; Y0; Z0; X0], options);
% plot results
uDose = zeros(1, numel(t));
for j = 1:numel(t)
[~, uDose(j)] = fM1_fun(t(j), y(j,:).', M1, i);
end
figure
plot(t, uDose)
xlabel('t')
ylabel('Amplitude')
title('Case 1: Approximation of $\delta(t - 10)$ with $\sigma = 0.001$', 'interpreter', 'latex', 'fontsize', 14)
figure
plot(t, y), grid on
legend('U', 'Y', 'Z', 'X')
xlabel('t')
ylabel('Amplitude')
title('Case 1: Responses of the Mice', 'interpreter', 'latex', 'fontsize', 14)
% Case 2
i = 2;
[t, y] = ode45(@(t, y) fM1_fun(t, y, M1, i), tspan, [U0; Y0; Z0; X0], options);
% plot results
uDose = zeros(1, numel(t));
for j = 1:numel(t)
[~, uDose(j)] = fM1_fun(t(j), y(j,:).', M1, i);
end
figure
plot(t, uDose)
xlabel('t')
ylabel('Amplitude')
title('Case 2: Approximation of $\delta(t - 10)$ with $\sigma = 0.005$', 'interpreter', 'latex', 'fontsize', 14)
figure
plot(t, y), grid on
legend('U', 'Y', 'Z', 'X')
xlabel('t')
ylabel('Amplitude')
title('Case 2: Responses of the Mice', 'interpreter', 'latex', 'fontsize', 14)
Hi @Júlia
Remark: The post was edited on Dec 10, 2025 to remove data and the main code per the OP's request.
Thank you for sharing the paper. However, I am unable to replicate the results. It is possible that my approximation of the Dirac Delta function is inaccurate. I only copied the data for the injected doses for [Subject 1] to test it in the code referenced in my previous comment. Please rest assured that I did not download your code.
The authors of the paper did not provide the numerical data, only displaying a representation of the data in an image (Figure 1). Moreover, I suspect that the authors used the discrete-time unit impulse sample 
instead of the approximation of the continuous-time Dirac Delta function 
.
instead of the approximation of the continuous-time Dirac Delta function .Note: This figure is not from the paper.
figure
plot(xData, yData)
xlabel('t, (Days)')
ylabel('Dose, (mg)')
title('Injected doses for Subject 1')
Answers (1)
Look at figure 3 (a) and (b) of the article. Each time the curve for U becomes discontinuous, a dosing is applied to the mice. If you compare this to figure 1, you will observe that the number of times dosing is applied is far greater than what figure (1) seems to suggest.
Thus what I would try if I were you is to take figure (3) and read out the time points when a jump appears and the height of the jump (which corresponds to Mi*ui in your code). This gives you two new arrays: "tdosei" and "jumpheighti". Then you can use these two arrays to repeat your computation.
If this procedure also gives results that far differ from the publication, you can additionally try to apply a simpler integration scheme with a fixed value of deltat = 0.2 days (although ode45 should of course be numerically more precise). I'm not sure whether the authors also used ode45. They report about a fixed value for the time step deltat = 0.2 days, but this could also mean that ode45 was forced to make outputs every 0.2 days (and computations were performed on a much finer grid). If the authors used a fixed time stepping scheme, this could also be a reason for results that differ from yours, and you could try such a simpler integration method.
Another thing that stroke me is that you use different arrays tDosei and tMassi. In my opinion, they must be equal because the weight of the mice was measured at the times when dosing was applied.
An indication whether you are on the right track is to check the tumor volume. It shouldn't increase far beyond 200 mm^3. In your present simulations, it reaches a maximum value of almost 1,2*10^6 mm^3 (for mouse 1) which is unrealistic.
I'd be interested if you could report about your future progress (hopefully).
1 Comment
when you suggest extracting the time points and jump heights directly from Figure 3, do you mean using a tool such as DataThief (as we did for Figure 1) to retrieve both the dosing times and the corresponding values?
I think if you don't want to contact the authors of the article, you don't have a better choice, do you ?
Categories
Find more on Code Generation in Help Center and File Exchange
on 10 Dec 2025
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
