6
$\begingroup$

Let me start by saying that this is my first post, so if there are suggestions on proper formatting that is annoying you, please let fly and I will do my best to modify.

I was doing well solving the topic in the title until I attempted to model second order reactions, which seemed to "break" Mathematica. I'll use an example of analogous code to get my point across.

So first I define my time series data as follows, including initial concentrations:

data = {{0., 182.115, 0.}, {1., 166.486, 1.87153}, {5., 136.178, 
   9.76618}, {10., 112.277, 19.1688}, {15., 111.448, 27.2978}, {25., 
   81.5446, 43.3609}, {35., 73.6962, 56.2542}, {45., 62.9892, 
   62.6874}, {55., 50.7161, 68.5274}, {65., 45.2139, 79.0054}};

{time, Xylose, Furfural} = Transpose[data];

Fin = Furfural[[1]];

Xin = Xylose[[1]];

Then I set up the various rate equations for the chemical species, and solve for the rate constants; so far, this is all well and good, and the output can be seen below.

sol = ParametricNDSolveValue[{X'[t] == -a X[t] - b X[t], X[0] == Xin, 
    F'[t] == a X[t] - c F[t], F[0] == Fin}, {X, F}, {t, 0, 70}, {a, b,
     c}];

transformeddata = {ConstantArray[Range@Length[{Xylose, Furfural}], 
      Length[time]] // Transpose, 
    ConstantArray[time, Length[{Xylose, Furfural}]], {Xylose, 
     Furfural}}~Flatten~{{2, 3}, {1}};

model[a_, b_, c_][i_, t_] := 
  Through[sol[a, b, c][t], List][[i]] /; 
   And @@ NumericQ /@ {a, b, c, i, t};

fit = NonlinearModelFit[transformeddata, 
   model[a, b, c][i, t], {a, b, c}, {i, t}];

fit["RSquared"]

fit["ParameterTable"]

Show[ListPlot[Table[Take[data, All, {1, l, l - 1}], {l, 2, 3}], 
  PlotStyle -> PointSize[0.02]], 
 Plot[{fit[1, t], fit[2, t]}, {t, 0, 70}]]

(This works great for modeling first order reaction kinetics)

The "Show stopper" comes when I try to represent a reaction with second order kinetics; all the above are first order kinetics. So for instance when I change the following:

sol = ParametricNDSolveValue[{X'[t] == **-a X[t]** - b X[t], X[0] == Xin, 
    F'[t] == **a X[t]** - c F[t], F[0] == Fin}, {X, F}, {t, 0, 70}, {a, b,
     c}];

into

sol = ParametricNDSolveValue[{X'[t] == **-a X[t] X[t]** - b X[t], X[0] == Xin, 
    F'[t] == **a X[t] X[t]** - c F[t], F[0] == Fin}, {X, F}, {t, 0, 70}, {a, b,
     c}];

then the output I get has the following errors:

Unsuccessful run with second order kinetics

In short, I'm not sure why this solve function is having trouble with representing second order reactions in this way. Any advice or information would be greatly appreciated.

$\endgroup$
3
  • $\begingroup$ Have a look at this question to see if the method provides any help to you. $\endgroup$ Oct 13, 2016 at 23:19
  • 3
    $\begingroup$ @bobthechemist meets @RobtheChemist?! $\endgroup$
    – Chris K
    Oct 14, 2016 at 0:53
  • 1
    $\begingroup$ @Chris, it was fated... $\endgroup$ Oct 14, 2016 at 2:28

1 Answer 1

5
$\begingroup$

You need to provide starting values for NonlinearModelFit.

First, I took your model and played with its parameters untill it looked good:

Show[Plot[model[0., 0.03, 0][1, t], {t, 0, 70}, PlotRange -> All], 
 ListPlot@Transpose@{time, Xylose}]

enter image description here

(You could wrap it in Manipulate.) It's very sensitive to the values of a, moderately sensitive to b and unnoticeably sensitive to c.

Next, I used those parameters as starting values (no errors at all):

fit = NonlinearModelFit[transformeddata, 
   model[a, b, c][i, t], {{a, 0.}, {b, 0.03}, {c, 0}}, {i, t}];

fit["RSquared"]

0.991887

fit["ParameterTable"]

enter image description here

Show[ListPlot[Table[Take[data, All, {1, l, l - 1}], {l, 2, 3}], 
  PlotStyle -> PointSize[0.02]], 
 Plot[{fit[1, t], fit[2, t]}, {t, 0, 70}]]

enter image description here

$\endgroup$
1
  • $\begingroup$ in this example, how can we define that the parameter should have positive value? $\endgroup$
    – martha
    Jan 12, 2021 at 11:05

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.