2
$\begingroup$

I try to fit parameters to ODE solver with chemical kinetic eqations:

R = 8.31; (* gas constant *)

(* inital contrations *)
c[v0_, T_, totalGasFlow_] := v0/totalGasFlow*101325/(R*T);
c0n2o = c[1.5, 298, 52];   
c0n2 = c[0, 298, 55];
c0o2 = c[0, 298, 52];

(* rate constants for reactions *)
calckarr[T_, nu_, eAct_] := 
 nu*Exp[-(eAct)/(R*T)]

kActFor[T_, a1f_] := calckarr[T, a1f, 111930];
kadsO2[T_, a3_] := calckarr[T, a3, 103246];
kdesO2[T_, a2_] := calckarr[T, a2, 117720];

(* rate equations of reactions*)
r1f[T_, a1f_] := kActFor[T, a1f]* cn2o[t]*(1 - theta[t]);
r2[T_, a2_] := kdesO2[T, a2]*(theta[t])^2;
r3[T_, a3_] := kadsO2[T, a3]*co2[t]*(1 - theta[t]);

(*rate of reactions for specific compounds *)
rtheta[T_, a1f_, a2_, a3_] := r1f[T, a1f] - 2 r2[T, a2] + 2 r3[T, a3];
rn2[T_, a1f_] := r1f[T, a1f]
rn2o[T_, a1f_] := -r1f[T, a1f]
ro2[T_, a2_, a3_] := r2[T, a2] - r3[T, a3]

(* ODE solver *)
sol[T_, a1f_, a2_, a3_] := NDSolve[
  {theta'[t] == rtheta[T, a1f, a2, a3],
   cn2'[t] == rn2[T, a1f],
   cn2o'[t] ==  rn2o[T, a1f],
   co2'[t] == ro2[T, a2, a3],
   theta[0] == 0,
   cn2[0] == c0n2,
   cn2o[0] == c0n2o,
   co2[0] == c0o2}, {cn2, cn2o, co2, theta}, {t, 10^-6, 3}]

(* experimental data *)
data={{374.15, 0.000627806}, {382.95, 0.00441}, {396.45, 0.00583}, {414.25,
   0.0061}, {433.35, 0.00821}, {451.95, 0.01279}, {470.15, 
  0.01131}, {488.05, 0.01178}, {507.15, 0.0122}, {526.65, 
  0.01351}, {546.35, 0.01393}, {566.05, 0.02557}, {586.45, 
  0.03422}, {606.45, 0.05179}, {626.75, 0.08718}, {647.05, 
  0.14891}, {664.75, 0.24017}, {682.75, 0.35589}, {700.05, 
  0.48896}, {717.15, 0.6169}, {734.55, 0.72361}, {750.95, 
  0.80726}, {765.55, 0.86888}, {781.45, 0.91729}, {797.25, 
  0.95264}, {812.55, 0.97429}}

I manipulated a1f, a2 and a3 to find first shots of values:

temps = Table[T, {T, 400, 800, 10}];
Manipulate[
 Show[ListPlot@data,
  ListLinePlot[
   Transpose@{temps, Table[
       (c0n2o - cn2o[0.4])/c0n2o /. sol[T, 10^a1f, 10^a2, 10^a3], {T, 
        temps}] // Flatten}, PlotStyle -> Red]],
 {{a1f,8.72}, 6, 14}, {{a2,10.9}, 6, 12}, {{a3,11.03}, 6, 12}]

enter image description here

So, the next step is fine-tuning and searching for perfect values with NonlinearModelFit and ParametricNDSolve:

conv = ParametricNDSolveValue[
  {theta'[t] == rtheta[T, a1f, a2, a3],
   cn2'[t] == rn2[T, a1f],
   cn2o'[t] ==  rn2o[T, a1f],
   co2'[t] == ro2[T, a2, a3],
   theta[0] == 0,
   cn2[0] == c0n2,
   cn2o[0] == c0n2o,
   co2[0] == c0o2},
 (c0n2o- cn2o[0.4])/c0n2o,
  {t, 10^-6, 3},
  {T, a1f, a2, a3}]

When I plot conv with parameters from Manipulate it works perfectly:

Show[Plot[conv[T, 10^8.72, 10^10.9, 10^11.03], {T, 400, 800}, 
  PlotStyle -> Red, PlotLegends -> {"model"}], 
 ListPlot[data, PlotLegends -> {"experimental data"}]]

enter image description here

But when I try to fit these values, even only for 1 parameter, it crashes:

fit = NonlinearModelFit[data, conv[T, a1f, 10^10.9, 10^11.03], a1f, T]


(* 
General::nomem: The current computation was aborted because there was insufficient memory available to complete the computation. 
OR
General::ovfl: Overflow occurred in computation. *)

Of course, memory is not the issue, I have 16GB of RAM and it fails after 2s of computations. Any idea what's going on?

Improve this question
$\endgroup$
3
  • $\begingroup$ When I run your code, the error that I get is InterpolatingFunction::dmval: Input value {0.4} lies outside the range of data in the interpolating function. Extrapolation will be used. In other words, it seems like NonlinearModelFit is searching outside the region of the interpolating function, and since interpolating functions are polynomials, they will blow up, sometimes quite rapidly, outside of the region on which the fit is made. I suspect that this is the problem. (My kernel also quits when I run your code.) A fix might be to add some constraints in NonlinearModelFit. Let me check. $\endgroup$
    – march
    Commented Sep 29, 2022 at 16:04
  • $\begingroup$ I can't seem to make anything work. I think there just might be an issue with trying to use a ParametricFunction, which amounts to an InterpolatingFunction once values have been entered, as the functional form in NonlinearModelFit. $\endgroup$
    – march
    Commented Sep 29, 2022 at 16:21
  • $\begingroup$ The extrapolation problem is known for me, it happens often with numerical calculations with differential equations. Still, NonlinearModelFit should reject values which blows the polynomials and start to search for the nearest place when it doesnt't crash. Maybe it is a problem which algorithm to use in such a situation. A simplier example of working fitting were solved by @MarcoB in this thread: link This question is just continuation. $\endgroup$
    – Lechuu
    Commented Sep 29, 2022 at 19:40

2 Answers 2

Reset to default
3
$\begingroup$
  • FindFit and Method -> NMinimize can find all of the {a1f, a2, a3}.
sol = FindFit[data, conv[T, a1f, a2, a3], {a1f, a2, a3}, T, 
  Method -> NMinimize]
Show[ListPlot[data], 
 Plot[conv[T, a1f, a2, a3] /. sol, {T, 0, 1000}, PlotStyle -> Red]]

enter image description here

  • Method -> "NMinimize" in NonlinearModelFit also work.
nlm = NonlinearModelFit[data, conv[T, a1f, a2, a3], {a1f, a2, a3}, T, 
  Method -> "NMinimize"]
Show[ListPlot[data], Plot[nlm[T], {T, 0, 1000}, PlotStyle -> Red]]
nlm["BestFitParameters"]

The same as FindFit.

enter image description here

enter image description here

Improve this answer
$\endgroup$
1
  • $\begingroup$ Wow, so simple, and works fantastic! Thank You for your help. $\endgroup$
    – Lechuu
    Commented Oct 4, 2022 at 16:05
2
$\begingroup$ 
$Version

(* "13.1.0 for Mac OS X x86 (64-bit) (June 16, 2022)" *)

Clear["Global`*"]

R = 831/100; (*gas constant*)

(*initial contrations*)
c[v0_, T_, totalGasFlow_] := v0/totalGasFlow*101325/(R*T);
c0n2o = c[3/2, 298, 52];
c0n2 = c[0, 298, 55];
c0o2 = c[0, 298, 52];

(*rate constants for reactions*)
calckarr[T_, nu_, eAct_] := nu*Exp[-(eAct)/(R*T)]

kActFor[T_, a1f_] := calckarr[T, a1f, 111930]
kadsO2[T_, a3_] := calckarr[T, a3, 103246]
kdesO2[T_, a2_] := calckarr[T, a2, 117720]

(*rate equations of reactions*)
r1f[T_, a1f_] := kActFor[T, a1f]*cn2o[t]*(1 - theta[t]);
r2[T_, a2_] := kdesO2[T, a2]*(theta[t])^2;
r3[T_, a3_] := kadsO2[T, a3]*co2[t]*(1 - theta[t]);

(*rate of reactions for specific compounds*)
rtheta[T_, a1f_, a2_, a3_] := 
  r1f[T, a1f] - 2 r2[T, a2] + 2 r3[T, a3];
rn2[T_, a1f_] := r1f[T, a1f]
rn2o[T_, a1f_] := -r1f[T, a1f]
ro2[T_, a2_, a3_] := r2[T, a2] - r3[T, a3]

(*experimental data*)
data = {
    {374.15, 0.000627806}, {382.95, 0.00441}, {396.45, 
     0.00583}, {414.25, 0.0061},
    {433.35, 0.00821}, {451.95, 0.01279},
    {470.15, 0.01131}, {488.05, 0.01178},
    {507.15, 0.0122}, {526.65, 0.01351},
    {546.35, 0.01393}, {566.05, 0.02557},
    {586.45, 0.03422}, {606.45, 0.05179},
    {626.75, 0.08718}, {647.05, 0.14891},
    {664.75, 0.24017}, {682.75, 0.35589},
    {700.05, 0.48896}, {717.15, 0.6169},
    {734.55, 0.72361}, {750.95, 0.80726},
    {765.55, 0.86888}, {781.45, 0.91729},
    {797.25, 0.95264}, {812.55, 0.97429}} //
   Rationalize[#, 0] &;

For ParametricNDSolve use

conv = ParametricNDSolve[{theta'[t] == rtheta[T, a1f, a2, a3], 
    cn2'[t] == rn2[T, a1f], cn2o'[t] == rn2o[T, a1f], 
    co2'[t] == ro2[T, a2, a3], theta[0] == 0, cn2[0] == c0n2, 
    cn2o[0] == c0n2o, co2[0] == c0o2}, {cn2, cn2o, co2, theta}, {t, 
    10^-6, 3}, {T, a1f, a2, a3}, WorkingPrecision -> 12];

(fit = NonlinearModelFit[data,
    (c0n2o - cn2o[T, a1f, 10^(109/10), 10^(1103/100)][
         2/5] /. conv)/c0n2o, a1f, T,
    WorkingPrecision -> 12])["BestFitParameters"]

(* NonlinearModelFit::sszero: The step size in the search has become less than the tolerance prescribed by the PrecisionGoal option, but the gradient is larger than the tolerance specified by the AccuracyGoal option. There is a possibility that the method has stalled at a point that is not a local minimum. *)

(* {a1f -> 5.54350617598*10^8} *)

Despite the warning,

Show[
 Plot[fit[T], {T, 400, 800},
  PlotStyle -> Red,
  PlotLegends -> {"model"}],
 ListPlot[data,
  PlotLegends -> {"experimental\ndata"}]]

enter image description here

Improve this answer
$\endgroup$
2
  • $\begingroup$ Thanks for Your help. For unknow reason, it doesn't work for me, maybe because of different version? I work on 12.1.1, Windows 10. The evaluation of the code takes more than minute for one variable. However, @cvgmt answer below works pretty well. $\endgroup$
    – Lechuu
    Commented Oct 4, 2022 at 16:05
  • $\begingroup$ Works with version 12.1.1 for Mac OS X x86 (64-bit) (June 19, 2020) The AbsoluteTiming for the plots was just under 22 seconds. $\endgroup$
    – Bob Hanlon
    Commented Oct 5, 2022 at 0:08

Your Answer

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

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