3
$\begingroup$

I call for your help with because I'm getting troubles to fit a model defined through a system of ODEs. The system of ODEs is as follows:

X'[t]:= m[t].X[t]
X[t_] := {h[t], r[t], rh1[t], rh2[t], rh3[t]}
m[t_] := {{-k1*r[t], 0, ki1, 0, 0}, {0, -k1*h[t], ki1, 0, 0}, {0,   k1*h[t], -(ki1 + k2), ki2, 0}, {0, 0, k2, -(ki2 + k3), ki3}

(It describes a chemical reaction consisting in three sequential steps, the first one is a bimolecular association and the following steps are unimolecular reactions). And the actual quantity I need to evaluate is:

F:= aF.X[t]
aF:={0,aR,aRH1,aRH2,aRH3}

(This quantity represents what I can measure experimentally). Well, following the directions found in the Wizard and in several forums I've tried the following.

Definition of the model:

Clear[modelo];
modelo[k1_?NumberQ, ki1_?NumberQ, k2_?NumberQ, ki2_?NumberQ, k3_?NumberQ, ki3_?NumberQ, aR_?NumberQ, aRH1_?NumberQ, aRH2_?NumberQ, aRH3_?NumberQ, Ht_?NumberQ, Rt_?NumberQ]:= (modelo[k1, ki1, k2, ki2, k3, ki3, aR, aRH1, aRH2, aRH3, Ht, Rt] = Module[{X, X0, m, vectorR, aF, sol, F}, First[{F,      {X[t_] := {h[t], r[t], rh1[t], rh2[t], rh3[t]};
   X0 := X[0] == {Ht, Rt, 0, 0, 0};       
   m[t_] := {{-k1*r[t], 0, ki1, 0, 0}, {0, -k1*h[t], ki1, 0, 
      0}, {0, k1*h[t], -(ki1 + k2), ki2, 0}, {0, 0, 
      k2, -(ki2 + k3), ki3}, {0, 0, 0, k3, -ki3}};
   aF := {0, aR, aRH1, aRH2, aRH3};
   sol := 
    NDSolve[{X'[t] == m[t].X[t], X0}, {h, r, rh1, rh2, rh3}, {t, 
      0, 500}];
   F := aF.Flatten[(X[t] /. sol)] - aR*Rt}}]])

The model so defined seems to work since, for example, it was evaluated and plotted smoothly by the following command:

aF = 0.0034;
Ht = 800;
Rt = 62; (*The last three quantities are not fitting parameters but fixed ones, whose values are known. *)
Fsim := With[{k1 = 0.0012, ki1 = 5, k2 = 0.43, ki2 = 0.0055, k3 = 0.01, ki3 = 0.0096, aRH1 = .032, aRH2 = .01, aRH3 = .012}, modelo[k1, ki1, k2, ki2, k3, ki3, aR, aRH1, aRH2, aRH3, Ht, Rt]];
LogLinearPlot[Evaluate[Fsim], {t, 0.001, 500}, PlotRange -> All, FrameLabel -> {"t(s)", "F"}]

Data and Fitting procedure:

To test that things work, I tried to fit the model to only one time series. Here is a fake data table:

data = Table[{t, .22 (1 - E^(-7.2 t)) + 0.10 (1 - E^(-0.084 t)) + 0.15 (1 - E^(-0.027 t))}, {t, 0, 500, 0.1}]

which actually proceeds from a fitting of that function to the real experimental data.

Finally, I've tried the following code to fit the model (only three parameters of it) to these data:

k1 = 0.0012;
k2 = 0.43;
ki2 = 0.0055;
k3 = 0.01;
ki3 = 0.0096;
aR = 0.0034;
Ht = 800;
Rt = 62; (*The last three quantities are not fitting parameters but fixed ones, whose values are known. *)
fit = NonlinearModelFit[data, {modelo[k1, ki1, k2, ki2, k3, ki3, aR, aRH1, aRH2, aRH3, Ht, Rt], {2 <= ki1 <= 20, 0.001 <= aRH1 <= 0.1, 0.001 <= aRH2 <= 0.1, 0.001 <= aRH3 <= 0.1}}, {{ki1, 5}, {aRH1, .032}, {aRH2, .01}, {aRH3, .012}}, t]

The problem is that it never worked, and I don't know where is the issue. It returns the following errors:

NonlinearModelFit::nrnum: The function value 1/2 ((-0.700176+<<3>>+0.0034 InterpolatingFunction[{{<<2>>}},{4,7,1,{<<1>>},{<<1>>},0,0,0,0,Automatic, },{},False},{{<<563>>}},{Developer`PackedArrayForm,{<<564>>},{<<1126>>}},{Automatic}][t])^2+(-0.700176+<<3>>+0.0034 InterpolatingFunction[{{<<2>>}},{4,7,1,{<<1>>},{<<1>>},0,0,0,0,Automatic,{},{},False},{{<<563>>}},<<1>>,{Automatic}][t])^2+(-0.700175+<<3>>+0.0034 <<1>>)^2+(<<1>>)^2+(<<1>>)^2+<<1>>^2+<<39>>+<<1>>+(<<1>>)^2+(<<1>>)^2+(<<1>>)^2+(-0.536879+<<3>>+0.0034 InterpolatingFunction[{{<<2>>}},<<3>>,{<<9>>}][t])^2+<<1>>) is not a real number at {ki1,aRH1,aRH2,aRH3} = {5.,0.032,0.01,0.012}. >>

All suggestions will be welcome.

$\endgroup$
  • $\begingroup$ Have you considered using ParametricNDSolve[] for this? $\endgroup$ – J. M. is away Aug 10 '15 at 5:34
  • $\begingroup$ Hi. No, I haven't. I thought the issue should not be in the NDsolve evaluation, because Evaluate[modelo[k1, ki1, k2, ki2, k3, ki3, aR, aRH1, aRH2, aRH3, Ht, Rt]], with a given set of parameter values, rendered the right solution. I will tested anyway. $\endgroup$ – jerefloyd Aug 10 '15 at 5:47
  • $\begingroup$ @Guess who it is. I've just tried your suggestion replacing the "sol" expression with sol := ParametricNDSolve[{X'[t] == m[t].X[t], X0}, {h, r, h1, rh2, rh3}, {t, 0, 500}, {k1, ki1, k2, ki2, k3, ki3, Ht, Rt}] and it didn't work. The evaluation of the model in this form failed when running Evaluate[With[{k1 = 0.0012, ki1 = 5, k2 = 0.43, ki2 = 0.0055, k3 = 0.01, ki3 = 0.0096, aRH1 = .032, aRH2 = .01, aRH3 = .012}, modelo[k1, ki1, k2, ki2, k3, ki3, aR, aRH1, aRH2, aRH3, Ht, Rt]]]. Thank you anyway. $\endgroup$ – jerefloyd Aug 10 '15 at 6:10
  • $\begingroup$ Welcome to Mathematica.SE! I suggest the following: 0) Browse the common pitfalls question 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Read the faq! 3) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! $\endgroup$ – Dr. belisarius Aug 10 '15 at 8:25
2
$\begingroup$

As I find your expressions difficult to work with, let me rewrite them somewhat:

Clear[modelo];
modelo[k1_?NumberQ, ki1_?NumberQ, k2_?NumberQ, ki2_?NumberQ, k3_?NumberQ, 
       ki3_?NumberQ, aR_?NumberQ, aRH1_?NumberQ, aRH2_?NumberQ, aRH3_?NumberQ, 
       Ht_?NumberQ, Rt_?NumberQ] := 
  modelo[k1, ki1, k2, ki2, k3, ki3, aR, aRH1, aRH2, aRH3, Ht, Rt] =
  Module[{X, m,  aF,  sol},
   X[t_] := {h[t], r[t], rh1[t], rh2[t], rh3[t]};
   m[t_] := {{-k1*r[t], 0, ki1, 0, 0}, {0, -k1*h[t], ki1, 0, 0}, 
             {0, k1*h[t], -(ki1 + k2), ki2, 0}, {0, 0, k2, -(ki2 + k3), ki3}, 
             {0, 0, 0, k3, -ki3}};
   aF = {0, aR, aRH1, aRH2, aRH3};
  sol = NDSolve[ {X'[t] == m[t].X[t], X[0] == {Ht, Rt, 0, 0, 0}}, X[t], {t, 0, 500}];
  Function[{tu}, Evaluate[aF.Flatten[X[t] /. sol] - aR*Rt] /. t :> tu]]

data = Table[{t, .22 (1 - E^(-7.2 t)) + 0.10 (1 - E^(-0.084 t)) + 
                  0.15 (1 - E^(-0.027 t))}, {t, 0, 500, 1}];

k1 = 0.0012; k2 = 0.43; ki2 = 0.0055; k3 = 0.01;  ki3 = 0.0096;
aR = 0.0034; Ht = 800; Rt = 62;
fit = 
 NonlinearModelFit[data, 
   {modelo[k1, ki1, k2, ki2, k3, ki3, aR, aRH1, aRH2, aRH3, Ht, Rt][t], 
    {2 <= ki1 <= 20, 0.001 <= aRH1 <= 0.1, 0.001 <= aRH2 <= 0.1, 0.001 <= aRH3 <= 0.1}}, 
   {{ki1, 5}, {aRH1, .032}, {aRH2, .01}, {aRH3, .012}}, t,  
  Method -> {NMinimize, Method -> "NelderMead"}]

fit["BestFitParameters"]
(* {ki1 -> 7.23261, aRH1 -> 0.0364504, aRH2 -> 0.0108148,  aRH3 -> 0.0116011} *)
$\endgroup$
  • $\begingroup$ I used "NelderMead" for speed. If you need another method to improve convergence or accuracy, that's a different matter. $\endgroup$ – Dr. belisarius Aug 10 '15 at 8:56
  • $\begingroup$ Thanks. I'm testing your suggestion right now. No errors so far, but -even with NelderMead and MaxIterations -> 10- it's taking a long time. Meanwhile I'll try to figure out the mean changes you made on the code. $\endgroup$ – jerefloyd Aug 10 '15 at 13:29
  • $\begingroup$ Thank you, finally the fitting worked out. But, as noticed before, its too slow. For instance, 109 seconds for a run with NelderMead and MaxIterations -> 2! That is a serious limitation for what I pretend to do, which -instead of 1 time series of 5000 points and fitting of 4 parameters- involves data consisting of ~ 8 time series of 14000 points each, and the simultaneous fitting of at least 7 parameters. Please let me know if you have any idea to speed up this task. $\endgroup$ – jerefloyd Aug 10 '15 at 19:08
  • 1
    $\begingroup$ @jerefloyd Ok,but these sites don't work that way. This isn't a forum. Here, you ask a question, get some answers, you vote on them and accept an answer (Please see my welcome comment about it under your question). Then if you have further questions you post another and so on. I know it sounds ridiculous at first, but this mechanics have its reasons, and they work. $\endgroup$ – Dr. belisarius Aug 10 '15 at 19:23
  • $\begingroup$ To be honest, I rather disagree. I found this page through Google, but - predictably - had the same problem with speed. If the intention of SE is to provide useful answers, then that intention has not been met. (And yes, I did find the related solution that did solve my problem.) $\endgroup$ – linkhyrule5 Mar 8 '17 at 3:33

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.