Group Abstract Group Abstract

Message Boards Message Boards

How to fit parameters in a system of coupled differential equations

Roc White
Posted 10 years ago

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:

Out[1]= 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:

Out[2]= FindMinimum[objfun[{k1, k2, k3}], {k1, k2, k3}]

What's the correct way to perform this fitting? Any suggestion would be great!
POSTED BY: Roc White
5 Replies
Sort By:
POSTED BY: Sam Carrettie
Roc White
Roc White
Posted 10 years ago

Thanks for your comments. Actually I've been quite aware that FindMinimum / NMinimize are purely numeric functions; what I didn't know is that I have to declare my function to take only numerical arguments for FindMinimum to work and that simply making sure that I don't pass in non-numeric arguments is not enough. Note that I haven't changed my data at all and they have always been just numbers. Intuitively, my assumption was that by being a numeric solver FindMinimum would only feed my function numbers, and furthermore that,

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

should not be necessary as the documentation for ParametricFunction specifies "A ParametricFunction object pfun is evaluated by using pfun[pvals] where pvals are explicit numerical values for the parameters pars." You would think since k1, k2, k3 only appear as pfun's arguments, they should be automatically treated as numeric.

In addition, note that in the example of my first reply, for a single parametric differential equation, I don't need to declare any arguments to be numeric.

Of course now seeing how one could make the code to work, I begin to piece together a better picture of what symbolic computation means in Mathematica, but a simple line in Details and Options about requirements on f would certainly have helped a lot.
POSTED BY: Roc White
Roc White
Roc White
Posted 10 years ago

I just came across this article: How do I use ?NumericQ to affect order of evaluation? which explains very well the cause of similar issues. Simply put, these problems would happen even if you feed only numerical arguments to your function. You actually need _?NumericQ to declare that the function takes only numerical arguments, thus changing the order of evaluation, in order to get rid of the problems. I believe this is not a matter of careless programming or RTFM type of thing; it's a missing warning of potential pitfalls in an otherwise very comprehensive documentation.
POSTED BY: Roc White
Roc White
Roc White
Posted 10 years ago
POSTED BY: Roc White
Roc White
Roc White
Posted 10 years ago

I just want to add to the discussion (without confusing the original post) that all the different routes 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 found the identical result:

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 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}}
POSTED BY: Roc White
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Tag limit exceeded
Note: Only the first five people you tag will receive an email notification; the other tagged names will appear as links to their profiles.
Reply Preview
Attachments
Remove
or Discard