Skip to main content

differential equations - Fails to apply NonlinearModelFit on a numerically evaluated model defined as a system of ODEs


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.



Answer



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} *)

Comments

Post a Comment

Popular posts from this blog

functions - Get leading series expansion term?

Given a function f[x] , I would like to have a function leadingSeries that returns just the leading term in the series around x=0 . For example: leadingSeries[(1/x + 2)/(4 + 1/x^2 + x)] x and leadingSeries[(1/x + 2 + (1 - 1/x^3)/4)/(4 + x)] -(1/(16 x^3)) Is there such a function in Mathematica? Or maybe one can implement it efficiently? EDIT I finally went with the following implementation, based on Carl Woll 's answer: lds[ex_,x_]:=( (ex/.x->(x+O[x]^2))/.SeriesData[U_,Z_,L_List,Mi_,Ma_,De_]:>SeriesData[U,Z,{L[[1]]},Mi,Mi+1,De]//Quiet//Normal) The advantage is, that this one also properly works with functions whose leading term is a constant: lds[Exp[x],x] 1 Answer Update 1 Updated to eliminate SeriesData and to not return additional terms Perhaps you could use: leadingSeries[expr_, x_] := Normal[expr /. x->(x+O[x]^2) /. a_List :> Take[a, 1]] Then for your examples: leadingSeries[(1/x + 2)/(4 + 1/x^2 + x), x] leadingSeries[Exp[x], x] leadingSeries[(1/x + 2 + (1 - 1/x...
Post a Comment
Read more

How to thread a list

I have data in format data = {{a1, a2}, {b1, b2}, {c1, c2}, {d1, d2}} Tableform: I want to thread it to : tdata = {{{a1, b1}, {a2, b2}}, {{a1, c1}, {a2, c2}}, {{a1, d1}, {a2, d2}}} Tableform: And I would like to do better then pseudofunction[n_] := Transpose[{data2[[1]], data2[[n]]}]; SetAttributes[pseudofunction, Listable]; Range[2, 4] // pseudofunction Here is my benchmark data, where data3 is normal sample of real data. data3 = Drop[ExcelWorkBook[[Column1 ;; Column4]], None, 1]; data2 = {a #, b #, c #, d #} & /@ Range[1, 10^5]; data = RandomReal[{0, 1}, {10^6, 4}]; Here is my benchmark code kptnw[list_] := Transpose[{Table[First@#, {Length@# - 1}], Rest@#}, {3, 1, 2}] &@list kptnw2[list_] := Transpose[{ConstantArray[First@#, Length@# - 1], Rest@#}, {3, 1, 2}] &@list OleksandrR[list_] := Flatten[Outer[List, List@First[list], Rest[list], 1], {{2}, {1, 4}}] paradox2[list_] := Partition[Riffle[list[[1]], #], 2] & /@ Drop[list, 1] RM[list_] := FoldList[Transpose[{First@li...
Post a Comment
Read more

front end - keyboard shortcut to invoke Insert new matrix

I frequently need to type in some matrices, and the menu command Insert > Table/Matrix > New... allows matrices with lines drawn between columns and rows, which is very helpful. I would like to make a keyboard shortcut for it, but cannot find the relevant frontend token command (4209405) for it. Since the FullForm[] and InputForm[] of matrices with lines drawn between rows and columns is the same as those without lines, it's hard to do this via 3rd party system-wide text expanders (e.g. autohotkey or atext on mac). How does one assign a keyboard shortcut for the menu item Insert > Table/Matrix > New... , preferably using only mathematica? Thanks! Answer In the MenuSetup.tr (for linux located in the $InstallationDirectory/SystemFiles/FrontEnd/TextResources/X/ directory), I changed the line MenuItem["&New...", "CreateGridBoxDialog"] to read MenuItem["&New...", "CreateGridBoxDialog", MenuKey["m", Modifiers-...
Post a Comment
Read more