0
$\begingroup$

I'm trying to estimate the parameters for the following pair of differential equations

meq = gamma*v[t]*m[t] - mu*m[t];
veq = v[t]*(k - epsilon - mu) - v[t]^2*k + m[t]*l - m[t]^2*l - 
   m[t]*v[t]*(k + l + gamma);
equations = {m'[t] == meq, v'[t] == veq};

where gamma, k, l and epsilon are parameters to be estimated, and mu is a known parameter with value 4.58*10^(-5). Also I have initial conditions given by

initials = {m[0] == 2.03, v[0] == 3.09};
mu = 4.58*10^(-5)  

and the data for the function m is the following

data = {{60, 2.13597}, {300, 2.27247}, {390, 2.29472}, {420, 2.40096}, {630, 
  2.59312}, {660, 2.84918}, {780, 2.93677}, {1020, 3.02945}, {1110, 
  3.04794}, {1140, 3.05796}, {1140, 3.08739}, {1380, 3.21218}, {1380, 
  3.2873}, {1500, 3.3347}, {1680, 3.44467}, {1710, 3.47574}, {1710, 
  3.48421}, {1830, 3.50433}}

I defined a parametric solution for the pair of equations with

pfun = ParametricNDSolveValue[Join[equations, initials], m, {t, 0, 2000}, 
{epsilon, k, l, gamma}]

And then I tried to use FindFit as

fit = FindFit[data, pfun[epsilon, k, l, gamma][t], {epsilon, k, l, gamma}, 
t]

However this method seems to be super inefficient as it seems to take unreasonable amount of time to complete the fit. Is there a way maybe to form the parametric function in a way that the FindFit works faster or use the FindFit more efficiently? Of course if you have another faster method in mind it's very welcome.

I also tried to solve this with some different programs than Mathematica and they seems to be working much faster but in those I have no idea what algorithms they use and also in them I have serious data issues.

$\endgroup$
  • 1
    $\begingroup$ "seems to be super inefficient" - note that you had not even attempted to provide starting guesses for epsilon, k, l, gamma; thus, FindFit[] just assigns them all a value of 1 as a guess and iterates accordingly, which you have seen to be intolerably slow. $\endgroup$ – J. M. will be back soon Oct 2 '18 at 11:05
  • $\begingroup$ @J.M.issomewhatokay. I also tried that but it seemed to be no help to start them from some values that should be closer to their "real" values. $\endgroup$ – Kplusn Oct 2 '18 at 11:11
2
$\begingroup$

Here you need to use a different model. In this model, too, not everything is smooth, but the result is clear and fast

meq = gamma*v[t]*m[t] - mu*m[t];
veq = v[t]*(k - epsilon - mu) - v[t]^2*k + m[t]*l - m[t]^2*l - 
   m[t]*v[t]*(k + l + gamma);
equations = {m'[t] == meq, v'[t] == veq};
initials = {m[0] == 2.03, v[0] == 3.09};
mu = 4.58*10^(-5);
data = {{60, 2.13597}, {300, 2.27247}, {390, 2.29472}, {420, 
    2.40096}, {630, 2.59312}, {660, 2.84918}, {780, 2.93677}, {1020, 
    3.02945}, {1110, 3.04794}, {1140, 3.05796}, {1140, 
    3.08739}, {1380, 3.21218}, {1380, 3.2873}, {1500, 3.3347}, {1680, 
    3.44467}, {1710, 3.47574}, {1710, 3.48421}, {1830, 3.50433}};


 model[epsilon_?NumberQ, k_?NumberQ, l_?NumberQ, 
   gamma_?NumberQ] := (model[epsilon, k, l, gamma] = 
    First[m /. 
      NDSolve[{v'[t] == 
         v[t]*(k - epsilon - mu) - v[t]^2*k + m[t]*l - m[t]^2*l - 
          m[t]*v[t]*(k + l + gamma), 
        m'[t] == gamma*v[t]*m[t] - mu*m[t], m[0] == 2.03, 
        v[0] == 3.09}, {m, v}, {t, 2000}]]);

 fit = 
 FindFit[data, 
  model[epsilon, k, l, gamma][
   t], {{epsilon, -4.8}, {k, 4.9}, {l, 0.}, {gamma, .1}}, t, 
  PrecisionGoal -> 4, AccuracyGoal -> 4]

(*Out[]= {epsilon -> -4.78195, k -> 4.89863, l -> -0.00336569, 
 gamma -> 0.0769293}*)

 Show[
 Plot[model[epsilon, k, l, gamma][t] /. fit, {t, 0, 2000}], 
 ListPlot[data, PlotStyle -> Orange]]

fig1

$\endgroup$
  • $\begingroup$ What do you mean not everything is smooth? And how is the model different? I see thar its fast but how is it different. $\endgroup$ – Kplusn Oct 3 '18 at 5:24
  • $\begingroup$ Not smooth because there is a message system that not everything is smooth. My code does not contain a parametric function that slows down the solution. $\endgroup$ – Alex Trounev Oct 3 '18 at 5:33
  • $\begingroup$ Okay, thank you. Still testing this trough more. $\endgroup$ – Kplusn Oct 3 '18 at 5:59
  • $\begingroup$ How did you choose the initial quesses {{epsilon, -4.8}, {k, 4.9}, {l, 0.}, {gamma, .1}} for the four parameters? Because the ones you selected seem to be much better than other ones that I'm trying. $\endgroup$ – Kplusn Oct 3 '18 at 13:41
  • $\begingroup$ @Kplusn I found them in two steps. First, I took all the parameters equal to 1, except for l = 0. Got a solution and took it as the next step. $\endgroup$ – Alex Trounev Oct 3 '18 at 14:18

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.