5
$\begingroup$

Problem description: I'm studying a chemical reaction that can be modeled by a system of coupled differential equations, and I want to use measured data to determine the parameters that appear in the equations.

(1) I come up with the following code:

data = {{0, 0.269036323, 0, 0},
    {1.855, 0.26559627, 0.001414574, 0.000317798},
    {2.715, 0.265004681, 0.002081772, 0.000435464},
    {4.004, 0.26092304, 0.003181524, 0.000689863}}\[Transpose];
ti = data[[1, All]]; (* independent variable *)
ci = data[[2 ;; 4, All]]; (* three dependent variables *)
pfun = ParametricNDSolveValue[{
   c1'[t] == -k1/(1 + k3*c1[t]) - 2*k2*k3*c1[t]/(1 + k3*c1[t]),
   c2'[t] == k2*k3*c1[t]/(1 + k3*c1[t]),
   c3'[t] == k1/(1 + k3*c1[t]),
   c1[ti[[1]]] == ci[[1, 1]],
   c2[ti[[1]]] == ci[[2, 1]],
   c3[ti[[1]]] == ci[[3, 1]]},
  {c1, c2, c3}, {t, 0, 20}, {k1, k2, k3}];
FindMinimum[Sum[Total[(ci[[i, All]] - Map[pfun[k1, k2, k3][[i]], ti])^2], {i, 1, 3}],
  {{k1, 3.*^-4}, {k2, 1.74*^-3}, {k3, 3.81}}]

and the result is:

FindMinimum[\!\(
\*UnderoverscriptBox[\(\[Sum]\), \(i = 1\), \(3\)]\(Total[
\*SuperscriptBox[\((ci[[i, All]] - \(pfun[k1, k2, k3]\)[[i]] /@ 
       ti)\), \(2\)]]\)\), {{k1, 0.0003}, {k2, 0.00174}, {k3, 3.81}}]

It merely rephrased what I entered without evaluating it.

(2) I tried also setting up the equations this way:

objfun[k_] := Module[{},
  fun = NDSolveValue[{
     c1'[t] == -k[[1]]/(1 + k[[3]]*c1[t]) - 2*k[[2]]*k[[3]]*c1[t]/(1 + k[[3]]*c1[t]),
     c2'[t] == k[[2]]*k[[3]]*c1[t]/(1 + k[[3]]*c1[t]),
     c3'[t] == k[[1]]/(1 + k[[3]]*c1[t]),
     c1[ti[[1]]] == ci[[1, 1]],
     c2[ti[[1]]] == ci[[2, 1]],
     c3[ti[[1]]] == ci[[3, 1]]},
    {c1, c2, c3}, {t, 0, 20}];
  Sum[Total[(ci[[i, All]] - Map[fun[[i]], ti])^2], {i, 1, 3}]
  ];
FindMinimum[objfun[{k1, k2, k3}], {k1, k2, k3}]

Still it didn't evaluate:

FindMinimum[objfun[{k1, k2, k3}], {k1, k2, k3}]

What's the correct way to perform this fitting? Any suggestion would be great!

$\endgroup$
  • $\begingroup$ It seems that the solutions in this post may work, but it will still help if someone can explain why my code doesn't. Thanks! $\endgroup$ – Peng Bai Apr 26 '15 at 1:08
2
$\begingroup$

Define this function:

f[k1_?NumericQ, k2_?NumericQ, k3_?NumericQ] := 
 Sum[Total[(ci[[i, All]] - Map[pfun[k1, k2, k3][[i]], ti])^2], {i, 1, 3}] // Quiet

Then,

fit = NMinimize[f[k1, k2, k3], {k1, k2, k3}];
params = fit // Last
(*{k1 -> 0.000194805, k2 -> 0.0291469, k3 -> 0.109229}*)

Plot it,

Table[Show[
  ListPlot[Transpose[{ti, ci[[i]]}]],
  Plot[(pfun[k1, k2, k3] /.params)[[i]][t], {t, ti[[1]], ti // Last}
   ]], {i, 1, 3}]

enter image description here

$\endgroup$
  • $\begingroup$ @Roc White : Mathematica is like a retarded kid with an axe.The retarded kid will try to destroy everything with the axe, no matter what. The axe is called symbolic processing. So, sometimes, you have to stop it from doing this and declaring Numeric parameters sometimes solve the situation, and takes away the axe from the retarded kid. $\endgroup$ – Ivan Apr 26 '15 at 16:51
  • $\begingroup$ It's actually related to the order of evaluation, as explained here. One wants to prevent Mathematica trying to evaluate f[k1,k2,k3] symbolically upon encountering it in FindMinimum or NMinimize, although FindMinimum only calls f with numerical arguments. $\endgroup$ – Peng Bai May 20 '15 at 20:03
4
$\begingroup$

I just wanted to point out that all the routes above will work if there is only one differential equation:

data = NDSolveValue[{
     x''[t] - k1*(1 - x[t]^2)*x'[t] + k2*x[t] == 0, 
     x[0] == 2, x'[0] == 0} /. {k1 -> 1., k2 -> 1.}, 
   Table[{t, x[t] + RandomReal[{-.3, .3}]}, {t, 0, 10, .2}], {t, 10}];
dataT = data\[Transpose];
ti = dataT[[1, All]];
di = dataT[[2, All]];
pfun = ParametricNDSolveValue[{
    x''[t] - k1*(1 - x[t]^2)*x'[t] + k2*x[t] == 0,
    x[0] == 2, x'[0] == 0},
    x, {t, 0, 10}, {k1, k2}];

FindFit finds the best-fit parameters sucessfully:

fit = FindFit[data, pfun[k1, k2][t], {{k1, 2}, {k2, 0}}, t]
Out[1]= {k1 -> 1.09028, k2 -> 1.02729}

FindMinimum finds them too:

FindMinimum[Total[(di - Map[pfun[k1, k2], ti])^2], {k1, k2}]
Out[2]= {1.41041, {k1 -> 1.09028, k2 -> 1.02729}}

And the Module approach also produced the same result:

objfun[k_] := Module[{},
  fun = NDSolveValue[{
     x''[t] - k[[1]]*(1 - x[t]^2)*x'[t] + k[[2]]*x[t] == 0,
     x[0] == 2, x'[0] == 0},
    x, {t, 0, 10}];
  Total[(di - Map[fun, ti])^2]
  ]
FindMinimum[objfun[{k1, k2}], {{k1, 1.0}, {k2, 0.0}}]
Out[3]= {1.41041, {k1 -> 1.09028, k2 -> 1.02729}}
$\endgroup$
1
$\begingroup$

Going through this solution, it looks like that function f in FindMinimum[f,x] should be declared to have numerical parameters. Ivan's answer indicated this too.

Essentially, for the parametric function route to work, one needs:

func[k1_, k2_, k3_] := Total[(ci - Through[pfun[k1, k2, k3][ti]])^2, 2]
   /; And @@ NumericQ /@ {k1, k2, k3};
FindMinimum[func[k1, k2, k3], {{k1, 3.*^-4}, {k2, 1.74*^-3}, {k3, 3.81}}]

This is just another way (besides what Ivan suggested) to declare k1, k2, k3 to be numerical, but I've gotten rid of Sum in favor of Through.

For the Module route to work, just add a similar condition:

objfun[k_] := Module[{},
fun = NDSolveValue[{
   c1'[t] == -k[[1]]/(1 + k[[3]]*c1[t]) - 
     2*k[[2]]*k[[3]]*c1[t]/(1 + k[[3]]*c1[t]),
   c2'[t] == k[[2]]*k[[3]]*c1[t]/(1 + k[[3]]*c1[t]),
   c3'[t] == k[[1]]/(1 + k[[3]]*c1[t]),
   c1[ti[[1]]] == ci[[1, 1]],
   c2[ti[[1]]] == ci[[2, 1]],
   c3[ti[[1]]] == ci[[3, 1]]},
  {c1, c2, c3}, {t, 0, 20}];
Total[(ci - Through[fun[ti]])^2, 2]
] /; And @@ NumericQ /@ k;
$\endgroup$

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.