0
|
7735 Views
|
4 Replies
|
2 Total Likes
View groups...
Share

# 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
4 Replies
Sort By:
Posted 11 years ago
 Great!I will try this approach. Thanks!
Posted 11 years ago
 @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 11 years ago
 No that is a typo. I've removed.Apologies.
Posted 11 years ago
 Did you mean to have adjacent commas in kparam[t_?NumberQ,,ka_?NumberQ, sig_]  ...?
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.