Group Abstract Group Abstract

Message Boards Message Boards

Global data fitting for a chemical kinetics problem?

Dear all The problem of global data analysis is a common problem among experimenters, especially in chemical kinetics. I have different kinetics (variation of absorbance as a function of time) at different wavelengths. I would like to adjust a set of 10 kinetics. These kinetics are governed by the same kinetic law, the same rate constants (here k3,k4,k5 ) but the absorption coefficients (ea, eh and ef) depend on the wavelength. I usually do these fits with Igor (Wavemetrics) but without any analytical solution I turned to Mathematica. I tried the MultiNonlinearModelFit function 

Clear[k3, k4, k5, ea1, eh1, ef1, ea2, eh2, ef2]
sol = ParametricNDSolve[{Derivative[1][ag][t] == -k3*ag[t], 
    Derivative[1][gox][t] == k3*ag[t], 
    Derivative[1][h][t] == 
     k3*ag[t] - 2*k4*h[t]*h[t] + 2*k5*ff[t]*ff[t], 
    Derivative[1][ff][t] == 2*k4*h[t]*h[t] - 2*k5*ff[t]*ff[t], 
    ag[0] == 0.000006, h[0] == 0, gox[0] == 0, ff[0] == 0}, {ag, gox, 
    h, ff}, {t, 0, 2000}, {k3, k4, k5}];
{agFunc, goxFunc, hFunc, fFunc} = {ag, gox, h, f} /. sol;
model1[k3_, k4_, k5_, ea1_, eh1_, ef1_][
   t_] := (((agFunc[k3, k4, k5][t])*ea1) + (eh1*
      hFunc[k3, k4, k5][t]) + (ef1*fFunc[k3, k4, k5][t]));
model2[k3_, k4_, k5_, ea2_, eh2_, ef2_][
   t_] := (((agFunc[k3, k4, k5][t])*ea2) + (eh2*
      hFunc[k3, k4, k5][t]) + (ef2*fFunc[k3, k4, k5][t]));
fit = MultiNonlinearModelFit[{dat1, 
    dat2}, {model1[k3, k4, k5, ea1, eh1, ef1][t], 
    model2[k3, k4, k5, ea2, eh2, ef2][t]}, {k3, k4, k5, ea1, eh1, ef1, 
    ea2, eh2, ef2}, {t}];
fit["BestFitParameters"]

I have a concern with the output I don't understand. 

ParametricNDSolve::ndsz: At t == 3.127085338387159`, step size is effectively zero; singularity or stiff system suspected.

InterpolatingFunction::dmval: Input value {3.2} lies outside the range of data in the interpolating function. Extrapolation will be used.

...

NonlinearModelFit::nrlnum: The function value {-0.039527+1. f[1.,1.,1.][0.2],-0.0396433+1. f[1.,1.,1.][0.4],-0.0396255+1. f[1.,1.,1.][0.6],-0.0394311+1. f[1.,1.,1.][0.8],-0.0393989+1. f[1.,1.,1.][1.],-0.0394361+1. f[1.,1.,1.][1.2],-0.0394022+1. f[1.,1.,1.][1.4],-0.039186+1. f[1.,1.,1.][1.6],-0.0392537+1. f[1.,1.,1.][1.8],-0.0392642+1. f[1.,1.,1.][2.],<<32>>,-0.0384155+1. f[1.,1.,1.][8.6],-0.0383621+1. f[1.,1.,1.][8.8],-0.0385297+1. f[1.,1.,1.][9.],-0.0384413+1. f[1.,1.,1.][9.2],-0.0384188+1. f[1.,1.,1.][9.4],-0.0382994+1. f[1.,1.,1.][9.6],-0.0383953+1. f[1.,1.,1.][9.8],-0.0381658+1. f[1.,1.,1.][10.],<<850>>} is not a list of real numbers with dimensions {900} at {k3,k4,k5,ea1,eh1,ef1,ea2,eh2,ef2} = {1.,1.,1.,1.,1.,1.,1.,1.,1.}.

Yet it provides me with parameters which for some are far from the expected values: {k3 -> 76620.6, k4 -> -23293.5, k5 -> 10237.6, ea1 -> 169.232, eh1 -> 81.862, ef1 -> 14461.1, ea2 -> 157.347, eh2 -> 80.0374, ef2 -> 14761.5}

Parts of the data are 

dat2 = {{0.2, 0.0366929}, {0.4, 0.0368321}, {0.6, 0.0368786}, {0.8, 
   0.0367551}, {1, 0.036677}, {1.2, 0.036687}, {1.4, 0.0365921}, {1.6,
    0.0365438}, {1.8, 0.0364348}, {2, 0.0364835}, {2.2, 
   0.036489}, {2.4, 0.0363618}, {2.6, 0.0361545}, {2.8, 
   0.0362492}, {3, 0.0362564}, {3.2, 0.0362539}, {3.4, 
   0.0361269}, {3.6, 0.0362652}, {3.8, 0.0361286}, {4, 
   0.0361022}, {4.2, 0.0360444}, {4.4, 0.0360916}, {4.6, 
   0.0360261}, {4.8, 0.0360464}, {5, 0.0359304}, {5.2, 
   0.0358969}, {5.4, 0.0359568}, {5.6, 0.036002}, {5.8, 
   0.0358563}, {6, 0.0360431}, {6.2, 0.0358403}, {6.4, 
   0.0358596}, {6.6, 0.0358963}, {6.8, 0.0359946}, {7, 
   0.0358273}, {7.2, 0.0359167}, {7.4, 0.0357477}, {7.6, 
   0.0357798}, {7.8, 0.0358253}, {8, 0.0358452}, {8.2, 
   0.0358645}, {8.4, 0.035836}, {8.6, 0.0358447}, {8.8, 
   0.0357238}, {9, 0.0358621}, {9.2, 0.0358738}, {9.4, 
   0.0357445}, {9.6, 0.0356904}, {9.8, 0.0357929}, {10, 0.0355983}}

and 

dat1 = {{0.2, 0.039533}, {0.4, 0.0396493}, {0.6, 0.0396315}, {0.8, 
   0.0394371}, {1, 0.0394049}, {1.2, 0.0394421}, {1.4, 
   0.0394082}, {1.6, 0.039192}, {1.8, 0.0392597}, {2, 
   0.0392702}, {2.2, 0.039207}, {2.4, 0.0391561}, {2.6, 
   0.0390186}, {2.8, 0.038951}, {3, 0.0390279}, {3.2, 
   0.0389623}, {3.4, 0.0389681}, {3.6, 0.0389333}, {3.8, 
   0.0387812}, {4, 0.0388705}, {4.2, 0.0388193}, {4.4, 
   0.0387575}, {4.6, 0.0386211}, {4.8, 0.0388169}, {5, 
   0.0387259}, {5.2, 0.038528}, {5.4, 0.0387643}, {5.6, 
   0.03884}, {5.8, 0.0385277}, {6, 0.0385109}, {6.2, 0.0386881}, {6.4,
    0.0385117}, {6.6, 0.038618}, {6.8, 0.0386693}, {7, 
   0.0385695}, {7.2, 0.0385362}, {7.4, 0.0384269}, {7.6, 
   0.038459}, {7.8, 0.038489}, {8, 0.0384817}, {8.2, 0.0384481}, {8.4,
    0.0384566}, {8.6, 0.0384215}, {8.8, 0.0383681}, {9, 
   0.0385357}, {9.2, 0.0384473}, {9.4, 0.0384248}, {9.6, 
   0.0383054}, {9.8, 0.0384013}, {10, 0.0381718}}

other question, but perhaps not far from the problem: is it possible to reduce the search area for the right parameters by specifying an existence interval for each parameter. E.g. 0.001<k3<0.01, 1<k4<400, 500<ea, etc?

Thank you in advance for your advice and help.

Hélène
POSTED BY: Hélène PASQUIER
5 Replies
Sort By:
Jim Baldwin
Jim Baldwin, Retired
Posted 5 years ago
POSTED BY: Jim Baldwin

Indeed, you're right, I only sent you the beginning of the data : the data goes up to 850 and not 10.

I managed to do again what you did. Thank you so much for that breakthrough.I guess {k3, 0.22} means that the k3 values start from 0.22. Does it ?

Hélène
POSTED BY: Hélène PASQUIER
POSTED BY: Hélène PASQUIER

Dear Jim The version I used is the one I found at the link you sent me, that is to say

Options[MultiNonlinearModelFit] = {AccuracyGoal -> Automatic, 
   ConfidenceLevel -> 19/20, EvaluationMonitor -> None, 
   Gradient -> Automatic, MaxIterations -> Automatic,
    Method -> Automatic, PrecisionGoal -> Automatic, 
   StepMonitor -> None, Tolerance -> Automatic, 
   VarianceEstimatorFunction -> Automatic,
    Weights -> Automatic, WorkingPrecision -> Automatic, 
   "DatasetIndexSymbol" -> \[FormalN]
   }; (* \[Equal] Join[Options[NonlinearModelFit], \
{"DatasetIndexSymbol" -> \[FormalN]}] *)

MultiNonlinearModelFit[datasets_, form_, fitParams_, 
   independents : Except[_List], opts : OptionsPattern[]] := 
      MultiNonlinearModelFit[datasets, form, 
   fitParams, {independents}, opts];

MultiNonlinearModelFit[datasets_, form : Except[{__Rule}, _List], 
   fitParams_, independents_, opts : OptionsPattern[]] := 
      MultiNonlinearModelFit[
   datasets, <|"Expressions" -> form, "Constraints" -> True|>, 
   fitParams, independents, opts];

MultiNonlinearModelFit[
       datasets : {__?(MatrixQ[#1, NumericQ] &)}, 
       KeyValuePattern[{"Expressions" -> expressions_List, 
     "Constraints" -> constraints_}], 
       fitParams_List, 
       independents_List,
       opts : OptionsPattern[]
   ] := Module[{
        fitfun, weights,
        numSets = Length[expressions], 
        augmentedData = Join @@ MapIndexed[
               Join[ConstantArray[N[#2], Length[#1]], #1, 2] &,
               datasets
           ], 
        indexSymbol = OptionValue["DatasetIndexSymbol"]
    },
       Condition[
            fitfun = With[{

       conditions = 
        Join @@ ({#1, expressions[[#1]]} & ) /@ Range[numSets]
               }, 
                  Switch @@ Prepend[conditions, Round[indexSymbol]]
              ]; 
            weights = Replace[
                  OptionValue[Weights],
                  {
                       (list_List)?(VectorQ[#1, NumericQ] & ) /; 
         Length[list] === numSets :> 

        Join @@ MapThread[ConstantArray, {list, Length /@ datasets}], 

       list : {__?(VectorQ[#1, NumericQ] & )} /; 
         Length /@ list === Length /@ datasets :>
                            Join @@ list, 

       "InverseLengthWeights" :> 
        Join @@ (ConstantArray[1./#1, #1] & ) /@ Length /@ datasets
                   }
              ]; 
            NonlinearModelFit[
                 augmentedData,

     If[TrueQ[constraints], fitfun, {fitfun, constraints}], 
                 fitParams,
                 Flatten[{indexSymbol, independents}],
                 Weights -> weights, 

     Sequence @@ FilterRules[{opts}, Options[NonlinearModelFit]]
             ]
            ,
            numSets === Length[datasets]
        ]
   ];

Thanks for the typo, I had forgotten to correct it in my message. Hélène
POSTED BY: Hélène PASQUIER
Jim Baldwin
Jim Baldwin, Retired
Posted 5 years ago

What version do you have that has a function called MultiNonlinearModelFit?

Found it: https://resources.wolframcloud.com/FunctionRepository/resources/MultiNonlinearModelFit.

And you do have a typo:

{agFunc, goxFunc, hFunc, fFunc} = {ag, gox, h, f} /. sol;

should be

{agFunc, goxFunc, hFunc, fFunc} = {ag, gox, h, ff} /. sol;
POSTED BY: Jim Baldwin
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Tag limit exceeded
Note: Only the first five people you tag will receive an email notification; the other tagged names will appear as links to their profiles.
Reply Preview
Attachments
Remove
or Discard