Help implementing 4th-order Runge–Kutta for ODEs with Dirac impulses (discontinuous system)

49 visualizaciones (últimos 30 días)
Hello
  13 comentarios
Sam Chak
Sam Chak el 3 de Dic. de 2025
Editada: Sam Chak el 10 de Dic. de 2025 a las 11:39
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 .
Note: This figure is not from the paper.
figure
plot(xData, yData)
xlabel('t, (Days)')
ylabel('Dose, (mg)')
title('Injected doses for Subject 1')

Iniciar sesión para comentar.

Respuestas (1)

Torsten
Torsten el 3 de Dic. de 2025
Editada: Torsten el 4 de Dic. de 2025
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 comentario
Torsten
Torsten el 4 de Dic. de 2025
Editada: Torsten el 4 de Dic. de 2025
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 ?

Iniciar sesión para comentar.

Categorías

Más información sobre Code Generation en Help Center y File Exchange.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by