Message Boards Message Boards

Changing differential equation during NDSolve evaluation

Posted 11 years ago
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
POSTED BY: Pat Mac
4 Replies
Did you mean to have adjacent commas in
kparam[t_?NumberQ,,ka_?NumberQ, sig_]  ...
?
POSTED BY: Bruce Miller
Posted 11 years ago
No that is a typo.
I've removed.

Apologies.
POSTED BY: Pat Mac
@Pat, WhenEvent may be helpful. For instance, if the following function reaches steady state, a stimulus  is added: 

Code: 
sol=NDSolve[{x'[t] == 10-x[t],x[0]==1,WhenEvent[x'[t]<0.5,x[t]->100+x[t]]},x[t],{t,0,10}]
Plot[Evaluate[x[t]/.(sol[[1]])],{t,0,10}]
POSTED BY: Shenghui Yang
Posted 11 years ago
Great!

I will try this approach.

Thanks!
POSTED BY: Pat Mac
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract