Message Boards Message Boards

GROUPS:

Use FindFit to get the rate constants for a chemical reaction?

Posted 3 months ago
519 Views
|
3 Replies
|
3 Total Likes
|

Hi,

As the title says, I'm trying to find reaction rate constants using FindFit. Unfortunately, I'm getting the error "FindFit::nrlnum: The function value ... is not a list of real numbers with \ dimensions {512} at {Amp,k1,k1r,k3,k3r,k19,k20} = \ {424.5993497,82.,3.10^6,8.810^7,1.5*10^9,114000.,250.}"

The code is in the attached file. Please suggest ways to fix this - I don't understand what I'm doing incorrectly.

3 Replies

If you try to evaluate your model for numerical parameters, you'll notice it doesn't evaluate at all. This is because you have a typo in the definition of your model: there is one comma too many in between two of the k3r and k19 parameters. Once you fix that, you'll also notice that the model returns a list instead of a number, which is another problem you should fix (probably with First @ NDSolve[...]).

Finally, the way you call Model[187/100000, Amp, k1, k1r, k3, k3r, k19, k20][t] in FindFit will not work, since Model[...] will evaluate to an expression involving t and you then glue another t at the end. Instead, Model[...] should be returning a Function that accepts t as an argument. You can do that by replacing:

(Amp (InitialConditions - 1/2 (\[Beta][t] + \[Gamma][t]) +  3/2 \[Delta][t])

with:

Function[t, (Amp (InitialConditions - 1/2 (\[Beta][t] + \[Gamma][t]) +  3/2 \[Delta][t])]

and changing the solution variables in NDSolve from {\[Beta][t], \[Gamma][t], \[Delta][t], \[Zeta][t], \[Eta]} to {\[Beta], \[Gamma], \[Delta], \[Zeta], \[Eta]}. Doing so will make NDSolve return functions instead of expressions in t.

A few other remarks: it's usually better to use NumericQ for these kinds of matching patterns instead of NumberQ, since NumberQ has a few quirks that are often undesirable. For example, NumberQ[Pi] is False. Also, for fitting a function, you may want ask yourself if your model function should be using memoization, since FindFit shouldn't be asking for the same value twice anyway.

The whole model will now look like:

ClearAll[Model]
Model[InitialConditions_?NumericQ, Amp_?NumericQ, k1_?NumericQ, 
   k1r_?NumericQ, k3_?NumericQ, k3r_?NumericQ, k19_?NumericQ, 
   k20_?NumericQ] := 
  Function[t, Amp (InitialConditions - 1/2 (\[Beta][t] + \[Gamma][t]) +  3/2 \[Delta][t])] /. First@NDSolve[{
      \[Beta]'[t] == -(k1r \[Gamma][t] +  k3 (InitialConditions - 1/2 (\[Beta][t] + \[Gamma][t]) + 
               3/2 \[Delta][t])) \[Beta][t] +   k1 (InitialConditions - 1/2 (\[Beta][t] + \[Gamma][t]) + 
           3/2 \[Delta][t]) + k3r \[Delta][t],
       \[Gamma]'[t] ==  k1 (InitialConditions - 1/2 (\[Beta][t] + \[Gamma][t]) + 
           3/2 \[Delta][t]) - (k1r \[Beta][t] +  k19 ( M - \[Zeta][t] - \[Eta][t]) + 
           k20 \[Zeta][t]) \[Gamma][t],   \[Delta]'[t] == -k3r \[Delta][t] + 
        k3 \[Beta][ t] (InitialConditions - 1/2 (\[Beta][t] + \[Gamma][t]) +  3/2 \[Delta][t]),

      \[Zeta]'[t] ==   k19 ( M - \[Zeta][t] - \[Eta][t]) \[Gamma][t] -  k20 \[Zeta][t] \[Gamma][t],
      \[Eta]'[t] == k20 \[Zeta][t] \[Gamma][t],

      \[Beta][0] == Sqrt[(k1 InitialConditions)/k1r], \[Gamma][0] ==    Sqrt[(k1 InitialConditions)/k1r],          \[Delta][0] ==  0, \[Epsilon][0] == M, \[Zeta][0] == 0, 
      \[Eta][0] ==  0}, 
  {\[Beta], \[Gamma], \[Delta], \[Zeta], \[Eta]}, {t, 0, runtime}];

Test that it works:

Model[187/100000, Amp, k1, k1r, k3, k3r, k19, k20] /.  Rule @@@ {{Amp, 424.5993497}, {k1, 82}, {k1r, 3000000}, {k3, 8.8*10^7}, {k3r, 1.5*10^9}, {k19, 114000}, {k20, 250}}
%[1.]

This should return a function and a number.

Posted 3 months ago

Thank you, those corrections appear to make it work (though it is slow, so yes memoization probably should be used).

This should be in the documentation for NDSolve: "...changing the solution variables in NDSolve from {[Beta][t], [Gamma][t], [Delta][t], [Zeta][t], [Eta]} to {[Beta], [Gamma], [Delta], [Zeta], [Eta]}. Doing so will make NDSolve return functions instead of expressions in t."

The documentation for NumberQ and NumericQ should state the difference between the two functions that you mentioned: "It's usually better to use NumericQ for these kinds of matching patterns instead of NumberQ, since NumberQ has a few quirks that are often undesirable. For example, NumberQ[Pi] is False."

Thank you again.

You're welcome. I just checked, and you're right: the way the function is evaluated by FindFit, it's better to use memoization. I expected it to first evaluate the parameters and then insert the values for t (in which case you'd need to do the NDSolve only once), but it seems like that's not how FindFit substitutes the symbols.

As for NDSolve and NumberQ: these behaviors are documented in the Details sections and the examples, but you have to know where to look.

Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract