How to fit parameters in a system of coupled differential equations

Comments

Problem 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}]

Update It seems that the solutions in this post may work, but why not the code above?

Update 2 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.

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;

Update 4 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}}

Comments

Comments powered by Disqus