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

mathematical optimization - Minimizing using indices, error: Part::pkspec1: The expression cannot be used as a part specification

I want to use Minimize where the variables to minimize are indices pointing into an array. Here a MWE that hopefully shows what my problem is. vars = u@# & /@ Range[3]; cons = Flatten@ { Table[(u[j] != #) & /@ vars[[j + 1 ;; -1]], {j, 1, 3 - 1}], 1 vec1 = {1, 2, 3}; vec2 = {1, 2, 3}; Minimize[{Total@((vec1[[#]] - vec2[[u[#]]])^2 & /@ Range[1, 3]), cons}, vars, Integers] The error I get: Part::pkspec1: The expression u[1] cannot be used as a part specification. >> Answer Ok, it seems that one can get around Mathematica trying to evaluate vec2[[u[1]]] too early by using the function Indexed[vec2,u[1]] . The working MWE would then look like the following: vars = u@# & /@ Range[3]; cons = Flatten@{ Table[(u[j] != #) & /@ vars[[j + 1 ;; -1]], {j, 1, 3 - 1}], 1 vec1 = {1, 2, 3}; vec2 = {1, 2, 3}; NMinimize[ {Total@((vec1[[#]] - Indexed[vec2, u[#]])^2 & /@ R...
Post a Comment
Read more

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

What is and isn't a valid variable specification for Manipulate?

I have an expression whose terms have arguments (representing subscripts), like this: myExpr = A[0] + V[1,T] I would like to put it inside a Manipulate to see its value as I move around the parameters. (The goal is eventually to plot it wrt one of the variables inside.) However, Mathematica complains when I set V[1,T] as a manipulated variable: Manipulate[Evaluate[myExpr], {A[0], 0, 1}, {V[1, T], 0, 1}] (*Manipulate::vsform: Manipulate argument {V[1,T],0,1} does not have the correct form for a variable specification. >> *) As a workaround, if I get rid of the symbol T inside the argument, it works fine: Manipulate[ Evaluate[myExpr /. T -> 15], {A[0], 0, 1}, {V[1, 15], 0, 1}] Why this behavior? Can anyone point me to the documentation that says what counts as a valid variable? And is there a way to get Manpiulate to accept an expression with a symbolic argument as a variable? Investigations I've done so far: I tried using variableQ from this answer , but it says V[1...
Post a Comment
Read more