# 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. Attachments:
3 Replies
Sort By:
Posted 3 months ago
 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.
 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.