0
$\begingroup$

I am a newcomer to Mathematica (use ver 9.0).

I just want to fit the experimental data to a system of odes' by using NonlinearModelFit. The data, contains time in the first column and product concentration in the second. The model contains 4 species, species p1[t] is the one I want to fit the data to.

data = {{18.283, 0.0003365}, {39.415, 0.0003892}, {60.547, 
    0.00045}, {81.679, 0.0005631}, {102.811, 0.0006446}, {123.943, 
    0.0006756}, {145.075, 0.0007655}, {166.207, 0.0008306}, {187.339, 
    0.0008814}, {208.471, 0.000896}, {229.603, 0.0009816}, {250.735, 
    0.0011191}, {271.867, 0.0010797}, {292.999, 0.0011473}, {314.131, 
    0.0012265}, {335.263, 0.0013187}, {356.395, 0.0012809}, {377.527, 
    0.0013523}, {398.659, 0.0013841}, {419.791, 0.0014507}, {440.923, 
    0.0014885}, {462.055, 0.0014898}, {483.187, 0.0015491}, {504.319, 
    0.0015838}, {525.451, 0.0015651}, {546.583, 0.0015903}, {567.715, 
    0.0016328}, {588.847, 0.0016423}, {609.979, 0.0016149}, {631.111, 
    0.0016508}, {652.243, 0.0016296}, {673.375, 0.0016465}, {694.507, 
    0.0016776}, {715.639, 0.0016668}, {736.771, 0.0016685}, {757.903, 
    0.0016967}, {779.035, 0.001717}, {800.167, 0.0016997}, {821.299, 
    0.001721}};
tmax = Max[data[[All, 1]]];
e0 == 0.0002;
k21 == 0.001;
kh == 0.0001;
model = ParametricNDSolveValue[{s'[
      t] == -16  s[t] (0.0002 - x1[t] - x2[t]) + k21 x1[t] - kh  s[t],
    x1'[t] == 16  s[t] (0.0002 - x1[t] - x2[t]) - (k21 + k2) x1[t],
    x2'[t] == k2 x1[t] - k3 x2[t],
    p1'[t] == k2 x1[t] + kh s[t],
    s[0] == 0.002, x1[0] == 0, x2[0] == 0, p1[0] == 0}, 
   p1, {t, 0, tmax}, {k2, k3}];
fit = NonlinearModelFit[data, 
  model[k2, k3][t], {{k2, 0.4}, {k3, 0.7}}, {t, 0, tmax}]
plotfit = Plot[fit[t], {t, 0, tmax}, PlotRange -> {0, 0.0021}];
plotdata = ListPlot[data, PlotStyle -> PointSize[0.01]];
Show[plotfit, plotdata]

When I run the code, it produces a lot of messages:

NDSolve::deqn: Equation or list of equations expected instead of True in the first argument
General::ivar: 0 is not a valid variable.
ParametricNDSolve::dsvar: 0.016777965285714284` cannot be used as a variable.
General::ivar: 0.016777965285714284` is not a valid variable.

Would anyone please help me to fix it.

Improve this question
$\endgroup$
4
  • 2
    $\begingroup$ You probably used = instead of == earlier. That leaves a definition lingering. Just quit the kernel (to purge the old definitions) and try again. $\endgroup$
    – Chris K
    Jun 20, 2019 at 16:12
  • $\begingroup$ The last argument in fit=NonlinearModelFit[...] should be t instead of ` {t, 0, tmax}`. $\endgroup$
    – Ulrich Neumann
    Jun 20, 2019 at 19:43
  • $\begingroup$ I fixed all the above, but errors still appear: ParametricNDSolveValue::ndnum: Encountered non-numerical value for a derivative at t$317 == 0.`. >> General::stop: Further output of ParametricNDSolveValue::ndnum will be suppressed during this calculation. >>NonlinearModelFit::nrlnum: The function value {<<1>>} is not a list of real numbers with dimensions {39} at {k2,k3} = {0.4,0.7}. >> $\endgroup$
    – Irrrr
    Jun 21, 2019 at 10:24
  • $\begingroup$ By fixing all of the above, do you mean that you replaced e0 == 0.0002; k21 == 0.001; kh == 0.0001; with e0 = 0.0002; k21 = 0.001; kh = 0.0001;? Secondly, your data looks like it intersects the y-axis at ~0.00025, but your initial condition is p1[0]==0. That could make it a difficult to fit if the model is constrained to go through 0. $\endgroup$
    – Tim Laska
    Jun 21, 2019 at 11:50

1 Answer 1

Reset to default
2
$\begingroup$

The following fixes your Equal versus Set error and sets up a Manipulate to help set initial guesses for parameters.

data = {{18.283, 0.0003365}, {39.415, 0.0003892}, {60.547, 
    0.00045}, {81.679, 0.0005631}, {102.811, 0.0006446}, {123.943, 
    0.0006756}, {145.075, 0.0007655}, {166.207, 0.0008306}, {187.339, 
    0.0008814}, {208.471, 0.000896}, {229.603, 0.0009816}, {250.735, 
    0.0011191}, {271.867, 0.0010797}, {292.999, 0.0011473}, {314.131, 
    0.0012265}, {335.263, 0.0013187}, {356.395, 0.0012809}, {377.527, 
    0.0013523}, {398.659, 0.0013841}, {419.791, 0.0014507}, {440.923, 
    0.0014885}, {462.055, 0.0014898}, {483.187, 0.0015491}, {504.319, 
    0.0015838}, {525.451, 0.0015651}, {546.583, 0.0015903}, {567.715, 
    0.0016328}, {588.847, 0.0016423}, {609.979, 0.0016149}, {631.111, 
    0.0016508}, {652.243, 0.0016296}, {673.375, 0.0016465}, {694.507, 
    0.0016776}, {715.639, 0.0016668}, {736.771, 0.0016685}, {757.903, 
    0.0016967}, {779.035, 0.001717}, {800.167, 0.0016997}, {821.299, 
    0.001721}};
lp = ListPlot[data, PlotLegends -> {"data-p1"}];
(* Define constants *)
tmax = Max[data[[All, 1]]];
e0 = 0.0002;
k21 = 0.001;
kh = 0.0001;
(* System of ODEs *)
system = {Derivative[1][p1][t] == 0.0001` s[t] + k2 x1[t], 
   Derivative[1][s][t] == -0.0001` s[t] + 0.001` x1[t] - 
     16 s[t] (0.0002` - x1[t] - x2[t]), 
   Derivative[1][x1][t] == -(0.001` + k2) x1[t] + 
     16 s[t] (0.0002` - x1[t] - x2[t]), 
   Derivative[1][x2][t] == k2 x1[t] - k3 x2[t], p1[0] == 0, 
   s[0] == 0.002`, x1[0] == 0, x2[0] == 0};
vbles = {p1, s, x1, x2};
(* Note position of p1 *)
modelpos = 1;
parms = {k2, k3};
pfun = ParametricNDSolveValue[system, vbles, {t, 0, tmax}, parms];
model[p__Real][t_] := #[[modelpos]] &@Through[pfun[p][t], List]
(* Build Manipulate sliders *)
Manipulate[global = {k2, k3};
  Show[
  lp,
  Plot[Evaluate[model[k2, k3][t]], {t, 0, tmax}, 
       PlotRange -> {0., 0.8}, PlotLegends -> {"fitted-p1"}],
   ImageSize -> Large], 
  {{k2, 0.4}, 0.004, 40., Appearance -> "Labeled"}, {{k3, 0.7}, 
  0.007, 70., Appearance -> "Labeled"}]
Dynamic@global
(* {32.15, 0.007} *)

Manipulate Image

Now, we will use FindFit to estimate the parameters using the initial guess defined by Dynamic above.

(* Initial Guess Based on Slider Values *)
initguess = MapThread[List, {parms, First@Dynamic@global}]
fit = Quiet@
  FindFit[data, model[k2, k3][t], initguess, t, 
   Method -> {NMinimize, Method -> {"DifferentialEvolution"}}]
Show[lp, Plot[Evaluate[model[k2, k3][t] /. fit], {t, 0, tmax}, 
      PlotRange -> {0., 0.8}, PlotLegends -> {"fitted-p1"}],
  ImageSize -> Large]
(* {{k2, 32.15}, {k3, 0.007}} *)
(* {k2 -> 2.68866, k3 -> 0.209279} *)

Real data fit

The fit is okay, but the model has an initial condition of $p1(t=0)==0$ whereas the data does not appear to intersect $(0,0)$. You may want to reconsider your model. To illustrate the point, we can generate data with test parameters $\{0.4,0.7\}$ as given in the OP and use the same bad initial guess from the slider to show that FindFit can indeed recover the test parameters. 

(* Test values for parameters *)
testvals = {0.4, 0.7}
initguess
{sol, samples} = First@NDSolve[
     (system /. Thread[parms -> testvals])~
      Join~{WhenEvent[Mod[t, tmax/50] == 0, 
        Sow[{t, p1[t], s[t], x1[t], x2[t]}]]},
     vbles, {t, 0., tmax}, MaxStepSize -> 1] // Reap;
samples = First@samples;
simdata = samples[[All, {1, modelpos + 1}]];
simlp = ListPlot[simdata, PlotLegends -> {"Sim Data - p1"}];
fit = Quiet@
  FindFit[samples[[All, {1, 2}]], model[k2, k3][t], initguess, t, 
   Method -> {NMinimize, Method -> {"DifferentialEvolution"}}]
Show[simlp, Plot[Evaluate[model[k2, k3][t] /. fit], {t, 0, tmax}, 
      PlotRange -> {0., 0.8}, PlotLegends -> {"fitted-p1"}],
  ImageSize -> Large]
(* {0.4, 0.7} *)
(* {{k2, 32.15}, {k3, 0.007}} *)
(* {k2 -> 0.398766, k3 -> 0.702116} *)

Fit on simulated data

Improve this answer
$\endgroup$

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.