# Estimating parameters for pair of differential equations efficiently

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 == 2.03, v == 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.

• "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. – J. M. will be back soon Oct 2 '18 at 11:05
• @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. – Kplusn Oct 2 '18 at 11:11

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 == 2.03, v == 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 == 2.03,
v == 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]] • What do you mean not everything is smooth? And how is the model different? I see thar its fast but how is it different. – Kplusn Oct 3 '18 at 5:24
• 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. – Alex Trounev Oct 3 '18 at 5:33
• Okay, thank you. Still testing this trough more. – Kplusn Oct 3 '18 at 5:59
• 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. – Kplusn Oct 3 '18 at 13:41
• @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. – Alex Trounev Oct 3 '18 at 14:18