1
$\begingroup$

I'm trying to solve kinetic differential equation depended on time with three parameters: preexponential factor a, reaction order n and temperature T. My experimental data depends on temperature, not time, so my approach is: solve ParametricNDSolve with few temperature values, then collect values for a certain time t for all teperatures and then fit result to experiment. When I assume that preexp factor and reaction order as numbers, it works well:

temps = Table[T, {T, 300, 600, 25}]; (*temperature range*)
Ea = 90000;  (*activation energy in J/mol*)
R = 8.31;  (*gas constant*)

sol = ParametricNDSolve[{c'[t] == -a*Exp[-Ea/(R*#)]*c[t]^n, 
     c[0] == 1}, c, {t, 0, 10}, {a, n}] & /@ temps;

c1 = c[10^9, 1.5] /. sol; (*assuming a=10^9 and n=1.5*)
c2 = c1[[#]][10] /. sol[[#]] & /@ Range[Length[temps]]; (*collect values for time=10 at all temperatures*)

conv = (c2[[1]] - c2[[#]])/c2[[1]] & /@Range[Length[c2]]; (*construct conversion values at all temperatures*)
data = {temps, conv} // Transpose;
noisedata = {temps,data[[#, 2]] + RandomReal[{-0.02, 0.02}] & /@Range[Length[conv]]} // Transpose;

My problem is: how to construct the code to fit parameters a and n to noisedata ? I've tried the same approach:

c3 = c[a, n] /. sol;
c4 = c3[[#]][10] /. sol[[#]] & /@ Range[Length[temps]];
conv2 = (c4[[1]] - c4[[#]])/c4[[1]] & /@ Range[Length[c4]];
data2 = {temps, conv2} // Transpose;
f = Interpolation[data2] (* interpolate to obrain continous function *)

Interpolation with parameters as before plot the correct answer:

Show[Plot[f[T] /. {a -> 10^9, n -> 1.5}, {T, 300, 600}], 
 ListPlot[data]]

enter image description here

But it fails when it comes to fit:

fit = NonlinearModelFit[noisedata, f[T], {{a, 10^9}, {n, 1.4}}, T]

(* NonlinearModelFit::fmgz: Encountered a gradient that is effectively zero. The result returned
 may not be a minimum; it may be a maximum or a saddle point.

or

NonlinearModelFit::nrlnum: The function value {0.0123183 +0. I,(...)} is not a list of real numbers with
 dimensions {13} at {a,n} = {1.*10^13,1.4}. *)

Which is not true, because applying these number to Interpolate gives real results. Aby other possibilities to solve that puzzle?

$\endgroup$

1 Answer 1

1
$\begingroup$

I think the main problem here is that you are trying to fit to an interpolating function, which no longer explicitly depends on your parameters, instead of to a model function containing those parameters.

Below I've refactored your code rather significantly because I think you were making your life harder than it needs to be :-) The main change is the use of ParametricNDSolveValue to directly calculate the conversion rate at $t=10$, rather than doing that by hand afterwards. This also provides a parametric function object that can be used as a model in your fitting.

Let's start by generating some exact solutions, as you did:

Ea = 90000; R = 8.31;
c10 = 
  ParametricNDSolveValue[
    {c'[t] == -a Exp[-Ea/(R T)]c[t]^n, c[0] == 1},
    (c[0] - c[10])/c[0], {t, 0, 15}, {a, n, T}
  ]

Plot[c10[10^9, 1.5, T], {T, 300, 600}]

plot of conversion at t=10 for different temperatures

Let's fabricate some noisy data now:

noisy=
  Table[
    {T, c10[10^9, 1.5, T] + RandomReal[{-0.03, 0.03}]},
    {T, 300, 600, 25}
  ];

ListPlot[noisy, Frame -> True]

scatter plot of noisy data

Now let's fit the noisy data to the model provided by c10[a, n, T], where $a$ and $E_a$ are parameters to be fit, and $T$ is the independent variable:

fit = NonlinearModelFit[noisy, c10[a, n, T], {{a, 5 10^7}, {n, 4}}, T];
fit["BestFitParameters"]

(* Out: {a -> 1.02589*^9, n -> 1.53829} *)

Above I've used relatively "bad" starting values to show that the fit is somewhat resilient as well, at least with this relatively low amount of noise. Finally let's compare the fitted model with the noisy data it was fitted to, to show that it is a good fit:

Plot[fit[T], {T, 300, 600}, Epilog -> Point[noisy], Frame -> True]

plot: line = fit, points = data. Shows good fit

$\endgroup$
1
  • $\begingroup$ wow! It is really simplier that I thought it should be :) Thank You very much for Your help! Now I'll try to apply this method to my much more complicated problem with many parameters to fit but it looks like it should work :) $\endgroup$
    – Lechuu
    Jul 2, 2022 at 17:20

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.