Hi,
I am using NDSolve to fit a set of differential equations to a biochemical model. In the example for this discussion I am providing a practice model of radioactive decay, however the actual model that I need to fit is much larger. I would like to run the NDSolve computation for a set amount of time before adding an additional stimulus to the simulation, (i.e. run to steady state and then add signal). However the way that I implement this is computationally very inefficient. I'm hoping that someone here could help me make this code more efficient.
The code below uses a chi^2 function to calculate the distance between the model predictions and the data given in dexp at times texp.
The function (kparam[ka_?NumberQ, kb_?NumberQ, sig_] := If[25 < t, ka*sig, ka] ) adds a stimulus to the model after 25 seconds of simulation time. Using this function in the manner that I have listed below is prohibitively expensive for this and larger models. Is there a better way to implement this?
Thank You,
Patrick
Clear["Global`*"]
texp = {0.2, 2.2, 4.0, 5.0, 6.0, 8.0, 11.0, 15.0, 18.0, 26.0, 33.0, 39.0, 45.0};
dexp = {35., 25., 22.1, 17.9, 16.8, 13.7, 12.4, 7.5, 4.9, 4.0, 2.4, 1.4, 1.1};
sig = 2.0;
kparam[t_?NumberQ,ka_?NumberQ, sig_] := If[t > 25, ka*sig, ka]
chi2[ka_?NumberQ, kb_?NumberQ] := Block[{sol, A, B, c}, sol =
NDSolve[
{ A'[t] == -kparam[t, ka, sig]*A[t],
B'[t] == kparam[t, ka, sig]*A[t] - kb*B[t],
c'[t] == kb*B[t],
A[0] == 35.,
B[0] == 0.,
c[0] == 0.},
{A, B, c},
{t, 0., 100.}][[1]];
Apply[Plus, (dexp - (A[t] /. sol /. t -> texp))^2]];
chi2[0.1, 0.1] // Timing
(* This is the comparison without the additional function for the signal *)
chi22[ka_?NumberQ, kb_?NumberQ] := Block[{sol, A, B, c}, sol =
NDSolve[
{ A'[t] == -ka*A[t],
B'[t] == ka*A[t] - kb*B[t],
c'[t] == kb*B[t],
A[0] == 35.,
B[0] == 0.,
c[0] == 0.},
{A, B, c},
{t, 0., 100.}][[1]];
Apply[Plus, (dexp - (A[t] /. sol /. t -> texp))^2]];
chi22[0.1, 0.1] // Timing
