0
$\begingroup$

I'm trying to fit some data with a system of differential equations :

data={{6., 60.1536}, {9., 57.9807}, {12., 60.9089}, {15., 59.4291}, {18., 61.3227}, {21., 61.8788}, {24., 67.2192}, {27., 66.2767}, {30., 68.705}, {33., 68.2107}, {36., 70.8731}, {39., 68.7269}, {42., 
 73.2306}, {45., 72.3068}, {48., 74.8006}, {51., 72.1975}, {54., 76.577}, {57., 75.5894}, {60., 76.342}, {63., 75.5134}, {66., 
 77.47}, {69., 76.6854}, {72., 78.7422}, {75., 78.6074}, {78., 81.0158}, {81., 82.8521}, {84., 85.1395}, {87., 85.211}, {90., 
84.5157}, {93., 83.622}, {96., 88.1703}, {99., 85.6195}, {102., 86.8345}, {105., 86.5568}, {108., 87.5942}, {111., 88.3053}, {114., 
 88.3475}, {117., 89.3993}, {120., 91.7091}, {123., 89.7268}, {126., 90.6704}, {129., 89.7999}, {132., 90.369}, {135., 88.7787}, {138., 
90.3022}, {141., 89.8267}, {144., 91.2241}, {147., 91.2859}, {150., 92.992}, {153., 91.0079}, {156., 93.0784}, {159., 90.8868}, {162., 
92.7426}, {165., 92.757}, {168., 94.4202}, {171., 92.2914}, {174., 
90.3876}, {177., 89.3376}, {180., 89.814}, {183., 88.9134}, {186., 
89.7058}, {189., 91.642}, {192., 90.3205}, {195., 87.8566}, {198., 87.6065},            {201., 87.0403}, {204., 87.3344}, {207., 87.2313}, {210., 
 87.3705}, {213., 86.9135}, {216., 87.2684}, {219., 87.2989}, {222., 
   85.4766}, {225., 85.3534}, {228., 86.535}, {231., 86.1929}, {234., 
86.089}, {237., 85.9466}, {240., 85.1389}, {243., 85.0242}, {246., 
 84.4313}, {249., 83.7604}, {252., 81.9419}, {255., 83.773}, {258., 
 82.7046}, {261., 84.7331}, {264., 86.0393}, {267., 84.7472}, {270., 
  79.1677}, {273., 80.9426}, {276., 79.9624}, {279., 75.5272}, {282., 
 79.3103}, {285., 80.8015}, {288., 81.3927}, {291., 80.1678}, {294., 
 80.268}, {297., 79.9067}, {300., 76.9766}, {303., 81.8132}, {306., 
 73.6449}, {309., 76.4059}, {312., 76.4056}, {315., 81.7311}, {318., 
 80.8468}, {321., 80.958}, {324., 86.9248}, {327., 78.3434}, {330., 
 74.8752}, {333., 78.0912}, {336., 81.5165}, {339., 72.7919}, {342., 
  74.2966}, {345., 79.2233}, {348., 81.9791}, {351., 74.3276}, {354., 
 85.1221}, {357., 78.6944}, {360., 75.8183}, {363., 75.6696}, {366., 
 75.9147}, {369., 76.3326}, {372., 80.0048}, {375., 79.8311}, {378., 
 79.0427}, {381., 81.8084}, {384., 73.5742}, {387., 84.2291}, {390., 
 84.9122}, {393., 82.6657}, {396., 78.2888}, {399., 90.0235}, {402., 
 83.3667}, {405., 81.7737}, {408., 81.19}, {411., 82.3131}, {414., 
  79.8072}, {417., 74.4822}, {420., 75.6291}, {423., 82.2655}, {426., 
  73.704}, {429., 81.4184}, {432., 72.1127}, {435., 74.7053}, {438., 
  79.4664}, {441., 86.4491}, {444., 79.5096}, {447., 77.1761}, {450., 
  83.082}, {453., 80.3418}, {456., 85.3873}, {459., 85.7409}, {462., 
  73.3735}, {465., 72.2276}, {468., 82.7752}, {471., 71.6917}, {474., 
  78.5233}, {477., 82.4042}, {480., 83.8073}, {483., 91.5845}, {486., 
  82.8906}, {489., 87.3935}, {492., 89.9856}, {495., 74.1819}, {498., 
  77.5752}, {501., 82.6796}, {504., 79.2659}, {507., 81.5865}, {510., 
  82.709}, {513., 88.4083}, {516., 81.7317}, {519., 76.2638}, {522., 
  86.2863}, {525., 93.2163}, {528., 82.6943}}

And now the equations :

Rc = 50

 model[R0_?NumberQ, c_?NumberQ, d_?NumberQ] :=  Module[{y, m, x}, 
 First[y /. NDSolve[{y'[x] == (0.20 - 0.005*(y[x] - Rc) - m[x]*y[x]), 
  m'[x] == c - d*m[x], y[0] == R0, m[0] == 0}, {y, m}, {x, 0, 
  530}, Method -> {"EquationSimplification" -> "Residual"}]]]

 nlm = NonlinearModelFit[data, {model[59, c, d][x], 0 < c, 0 < d }, {{c, 0.5}, {d, 1}}, x]

But I've been waiting for a very long time and Mathematica didn't give me the solutions.

So I wanted to ask you if there was some "Method-> " I could add in my constrained NonlinearModelFit in order to get faster the result. I looked for such a method but what I read left me very unclear on how to do it exactly... Apparently most of the methods are not ok for constrained problems, and I tried Method->"NMinimize" which did not change anything. But apparently it has to be written : Method-> {"NMinimize",Method->...} and I don't know what to write exactly....

Thank you in advance

$\endgroup$
  • 1
    $\begingroup$ I think you need a better model or better starting values. If you plot the predictions from your current starting values (ListPlot[Transpose[{data[[All, 1]], model[59, 0.5, 1][#] & /@ data[[All, 1]]}]]), you won't get a curve that looks anything like the data. Also, for down the road, the variance about the curve is not constant (one of the assumptions with the defaults of NonlinearModelFit) but that can be examined after a reasonable fit (assuming constant variance) is found. $\endgroup$ – JimB Apr 15 '19 at 15:38
  • $\begingroup$ @JimB Thank you for those explanations. Could you explain me please what : ListPlot[Transpose[{data[[All, 1]], model[59, 0.5, 1][#] & /@ data[[All, 1]]}]] exactly means ? $\endgroup$ – J.A Apr 15 '19 at 15:46
  • 1
    $\begingroup$ The predictions for your starting values are found with model[59, 0.5, 1][#] & /@ data[[All, 1]]. The # symbol is replaced with each value of the predictor (data[[All,1]]). Then Transpose[{x,y}] constructs the predicted dataset. $\endgroup$ – JimB Apr 15 '19 at 16:01
  • 1
    $\begingroup$ Note that @CarlWoll 's answer suggests the current model is not a very good fit. In addition, looking at fit["CorrelationMatrix"] gives a result of $\left( \begin{array}{cc} 1. & 1. \\ 1. & 1. \\ \end{array} \right)$ which can suggest an overparameterized model or that lots of combinations of $c$ and $d$ give the same (poor) fit. $\endgroup$ – JimB Apr 15 '19 at 16:11
1
$\begingroup$

I would use ParametricNDSolveValue to create the model, because then parameter derivatives of the model can be used when creating the fit:

Rc = 50;
model = ParametricNDSolveValue[
    {
    y'[x] == 0.20 - 0.005 (y[x] - Rc) - m[x] y[x],
    m'[x] == c - d m[x],
    y[0] == R0,
    m[0] == 0
    },
    y,
    {x, 0, 530},
    {R0, c, d}
];

Then, use NonlinearModelFit:

fit = NonlinearModelFit[
    data,
    {model[59, c, d][x], 0 < c, 0 < d},
    {{c, .5}, {d, 1}},
    x
]; //AbsoluteTiming

{1.17995, Null}

Finally, plot the data and the fit:

Show[
    ListPlot[data],
    Plot[fit[x], {x, 0, 550}]
]

enter image description here

$\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.