Group Abstract Group Abstract

Message Boards Message Boards

How to fit parameters in a system of coupled differential equations

Roc White
Posted 11 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:

Thanks for actually coming back and posting the answer you found.

How could a new user such as me figure out this sort of stuff myself by using Mathematica's documentation?

If you mean that FindMinimum is purely numerical function, then going carefully through all the docs examples will give you a good idea, especially the Details and Properties & Relations sections. Once you get a sense that function is purely numeric then making sure you do not feed it symbolic values is matter of careful debugging. Note some functions can be both numeric and symbolic:

enter image description here

Not sure if this helps.
POSTED BY: Sam Carrettie
Roc White
Roc White
Posted 11 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 11 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 11 years ago

Found the solution, see here and here.

Basically, it looks like that function f in FindMinimum[f,x] should be declared to have numerical parameters. Therefore 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}}]

For the Module route to work, one needs to 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;

I'm still puzzled by how internal functions such as FindMinimum work. How could a new user such as me figure out this sort of stuff myself by using Mathematica's documentation? Thanks!
POSTED BY: Roc White
Roc White
Roc White
Posted 11 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