5
$\begingroup$

I need to fit a system of ODEs to data. The solutions which could work for me are given in the following posts

How to fit 3 data sets to a model of 4 differential equations? ODE fitting to dataset

However, when I reran both solutions proposed in those threads and apparently accepted by the community as working and valid, I got a series of errors (print after calling NonlinearModelFit). Could you please explain why? My understanding is that the model cannot return the values for A, B, or X depending on the value of the dummy variable i, but I don't know how to resolve the issue and why it appears in the first place.

The code follows below.

    a'[t] == -k1 a[t] b[t] + k2 x[t], a[0] == 1,
    b'[t] == -k1 a[t] b[t] + k2 x[t] - k3 b[t] x[t], b[0] == 1,
    x'[t] == k1 a[t] b[t] - k2 x[t] - k3 b[t] x[t], x[0] == 0
    }, {a, b, x}, {t, 0, 10}, {k1, k2, k3}
   ];

abscissae = Range[0., 10., 0.1];
ordinates = With[{k1 = 0.85, k2 = 0.15, k3 = 0.50},
   Through[sol[k1, k2, k3][abscissae], List]
  ];

data = ordinates + RandomVariate[NormalDistribution[0, 0.1^2], Dimensions[ordinates]];
ListLinePlot[data, DataRange -> {0, 10}, PlotRange -> All, AxesOrigin -> {0, 0}]

transformedData = {
    ConstantArray[Range@Length[ordinates], Length[abscissae]] // Transpose,
    ConstantArray[abscissae, Length[ordinates]],
    data
   } ~Flatten~ {{2, 3}, {1}};

model[k1_, k2_, k3_][i_, t_] := 
  Through[sol[k1, k2, k3][t], List][[i]] /;
    And @@ NumericQ /@ {k1, k2, k3, i, t};

fit = NonlinearModelFit[
   transformedData,
   model[k1, k2, k3][i, t],
   {k1, k2, k3}, {i, t}
  ];

ParametricNDSolveValue::ndsz: At t$2274999 == 2.683730722253005`, step size is effectively zero; singularity or stiff system suspected. >>
InterpolatingFunction::dmval: Input value {2.7} lies outside the range of data in the interpolating function. Extrapolation will be used. >>
$\endgroup$
5
  • 1
    $\begingroup$ Post you refer is dated as 2013, and in that time there been version 8 or 9. In every version some Automatic method used in NDSolve. Now in v.12 and 12.1 it is quite different method then in v. 8 and 9. What version you are running? $\endgroup$ Apr 18 '20 at 22:38
  • $\begingroup$ @AlexTrounev I am using 10.0.2.0 on Mac. Just to be sure, in your answer below you only added option Method in NonlinearModelFit? It now runs fine. $\endgroup$ Apr 18 '20 at 23:16
  • $\begingroup$ Yes, it runs fine in v.12.0.0. We can try several options but this one is good. $\endgroup$ Apr 18 '20 at 23:40
  • $\begingroup$ @AlexTrounev Regarding the method option in NonlinearModelFit, which ones are most suitable for problems with constraints? Is there any one which should be preferred (apart from Automatic)? $\endgroup$ Apr 19 '20 at 11:01
  • $\begingroup$ The method depends on constraints and equations we solve. It could be better you show us your particularly problem and we look how it could be solve. $\endgroup$ Apr 19 '20 at 11:14
6
$\begingroup$

We can use some options to solve this problem

sol = ParametricNDSolveValue[{a'[t] == -k1 a[t] b[t] + k2 x[t], 
    a[0] == 1, b'[t] == -k1 a[t] b[t] + k2 x[t] - k3 b[t] x[t], 
    b[0] == 1, x'[t] == k1 a[t] b[t] - k2 x[t] - k3 b[t] x[t], 
    x[0] == 0}, {a, b, x}, {t, 0, 10}, {k1, k2, k3}
   ];

abscissae = Range[0., 10., 0.1];
ordinates = 
  With[{k1 = 0.85, k2 = 0.15, k3 = 0.50}, 
   Through[sol[k1, k2, k3][abscissae], List]];

data = ordinates + 
   RandomVariate[NormalDistribution[0, 0.1^2], Dimensions[ordinates]];
ListLinePlot[data, DataRange -> {0, 10}, PlotRange -> All, 
 AxesOrigin -> {0, 0}]

transformedData = {ConstantArray[Range@Length[ordinates], 
      Length[abscissae]] // Transpose, 
    ConstantArray[abscissae, Length[ordinates]], data}~
   Flatten~{{2, 3}, {1}};

model[k1_, k2_, k3_][i_, t_] := 
  Through[sol[k1, k2, k3][t], List][[i]] /; 
   And @@ NumericQ /@ {k1, k2, k3, i, t};

fit = NonlinearModelFit[transformedData, 
   model[k1, k2, k3][i, t], {k1, k2, k3}, {i, t}, 
   Method -> "Gradient"];

Now we check and plot

Show[
 Plot[Evaluate[Table[fit[i, t], {i, 3}]], {t, 0, 10}, 
  PlotLegends -> {a, b, x}], 
 ListPlot[data, DataRange -> {0, 10}, PlotRange -> All, 
  AxesOrigin -> {0, 0}]]

Figure 1

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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