# How to fit parameters in a system of coupled differential equations

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!

• 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! Apr 26, 2015 at 1:08

## 3 Answers

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


• @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.
– Ivan
Apr 26, 2015 at 16:51
• 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. May 20, 2015 at 20:03

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


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;