6
$\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!

Improve this question
$\endgroup$
1
  • 1
    $\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$
    – P. B.
    Apr 26, 2015 at 1:08

3 Answers 3

Reset to default
3
$\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

Improve this answer
$\endgroup$
2
  • $\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, 2015 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$
    – P. B.
    May 20, 2015 at 20:03
6
$\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}}
Improve this answer
$\endgroup$
2
$\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;
Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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