0
$\begingroup$

my project right now consists of doing a data fit of this system of differential equations with 2 variables. Below is what I have so far but I keep running into problems with FindFit(keeps giving me some warnings such as NDSolve::mxst: Maximum number of 200 steps reached at the point t == 7.642417604774492`. What I have so far:

deltad = 0.35; deltav = 2.3; deltaz = 0.35; rho = 1.5;
alpha = 1;
r = 0.385;
s = 0.626;
uv = 5*10^10; \[Alpha] = 1; a = 1; gammay = 1;
ud = 1*10^6;
n = 3500;
sigmay = 0.00001;
sigmax = 0.00001;
fg = Piecewise[{{1*10^6, 0 <= t <= 2}}];
fg1 = Piecewise[{{5*10^10, 3 <= t <= 5}}];

data = {{0.000745712`, 99.23612972`}, {2.002982849`, 
205.1058697`}, {4.005219985`, 205.1058697`}, {6.007457122`, 
222.7508264`}, {8.009694258`, 169.8159564`}, {10.01193139`, 
134.5260431`}, {12.01416853`, 55.12373805`}, {14.06089983`, 
19.83382471`}, {16.0186428`, 2.188868043`}, {17.97638578`, 
2.188868043`}, {19.97862292`, 2.188868043`}, {21.98086005`, 
28.65630305`}, {24.02759135`, 46.30125972`}, {25.98533433`, 
134.5260431`}, {28.07655978`, 258.0407397`}};
first = r*x[t] - \[Kappa]*((x[t]*z[t])/(
x[t] + y[t] + z[t] + d[t])) - (\[Beta]*v[t]*x[t])/(
x[t] + y[t] + z[t] + d[t]);
second = s*y[t] + (\[Beta]*v[t]*x[t])/(x[t] + y[t] + z[t] + d[t]) - 
alpha*y[t] - \[Kappa]*(y[t]*z[t])/(x[t] + y[t] + z[t] + d[t]);
third = n*alpha*y[t] - deltav*v[t] + fg1;
fourth = sigmay*y[t] - deltad*d[t] + fg;
fifth = rho*d[t] - deltaz*z[t];
model[\[Kappa]_?NumberQ, \[Beta]_?
NumberQ] := (model[\[Kappa], \[Beta]] = 
First[y /. 
 NDSolve[{x'[t] == first, y'[t] == second, v'[t] == third, 
   d'[t] == fourth, z'[t] == fifth, x[0] == 1.5*10^8, y[0] == 0, 
   v[0] == 0, d[0] == 0, z[0] == 0}, {x, y, v, d, z}, {t, 0, 28}, 
  MaxSteps -> 200]])
fittedh = 
FindFit[data, {model[\[Kappa], \[Beta]][
t], \[Kappa] < \[Beta]}, {{\[Beta], 0.0001}, {\[Kappa], 0.001}}, 
t]

Unfortunately my output does not match the data. I would appreciate some help/ suggestions to make better.

Just in case, someone is interested, this is what my output was like

\[Kappa] = 0.001001524973637588`;
\[Beta] = 0.00010044137152890746`;
hj = NDSolve[{x'[t] == first, y'[t] == second, v'[t] == third,  d'[t] 
== fourth, z'[t] == fifth, x[0] == 1.5*10^8, y[0] == 0, v[0] == 0,
d[0] == 0, z[0] == 0}, {x, y, v, d, z}, {t, 0, 28.07655978`}]
ll = Show[ListPlot[data], 
Plot[Evaluate[((x[t] + y[t])/10000) /. {hj}], {t, 0, 28}]]

Failed attempt

$\endgroup$
  • $\begingroup$ You have left off the definition of deltav. $\endgroup$ – JimB Jul 10 '17 at 16:17
  • $\begingroup$ sorry about that @JimBaldwin. I just edited it. $\endgroup$ – peace96 Jul 10 '17 at 16:20
  • $\begingroup$ model[\[Kappa]_?NumberQ, \[Beta]_?NumberQ] := (zed=y /. NDSolve[{x'[t]==first, y'[t]==second, v'[t]==third, d'[t]==fourth, z'[t]==fifth, x[0]==1.5*10^8, y[0]==0, v[0]==0, d[0]==0, z[0]==0}, {x,y,v,d,z}, {t,0,28.1}, MaxSteps->400][[1]]; Print[Show[ListPlot[data], Plot[zed[t], {t,0,28.1}]]];Total[Map[Norm[#[[2]]-zed[#[[1]]]] &, data]]); NMinimize[{zzed=model[\[Kappa], \[Beta]], 0<\[Kappa]< .002 && 0<\[Beta]< .0002}, {\[Kappa], \[Beta]}, StepMonitor:>Print[{zzed,\[Kappa],\[Beta]}]] shows your starting estimates are bad or your model is bad because the error never gets below 1000. $\endgroup$ – Bill Jul 10 '17 at 20:14
  • $\begingroup$ model[\[Kappa]_?NumberQ, \[Beta]_?NumberQ] := (zed = x /. NDSolve[{x'[t] == first, y'[t] == second, v'[t] == third, d'[t] == fourth, z'[t] == fifth, x[0] == 1500, y[0] == 0,v[0] == 0, d[0] == 0, z[0] == 0}, {x, y, v, d, z}, {t, 0, 28.1}, MaxSteps -> 400][[1]];Print[Show[ListPlot[data], Plot[(zed[t])/10, {t,0,28.1}]]];Total[Map[Norm[#[[2]] - (zed[#[[1]]] + zed[[[2]]])/10]&,data]]);NMinimize[{zzed = model[\[Kappa], \[Beta]], 0 < \[Kappa] < 0.0002 && 0 < \[Beta] < .005 }, [Beta]},StepMonitor :> Print[{zzed, \[Kappa], \[Beta]}]]. Is there an obvious reason why the fit always falls on t=5? $\endgroup$ – peace96 Jul 11 '17 at 2:18
  • $\begingroup$ Your data is {{T1,Y1},{T2,Y2}...}, the solution of your de is y[t]. For each {Ti,Yi} I calculated Norm[Yi-y[Ti]]. See how this was adding up the difference between your data points and the "solution" and that difference is what you want to minimize? You changed this to Norm[Yi-(y[Ti]+zed[[[2]]])/10] and I can't even guess what zed[[[2]]] is. So I'm sorry, but I have no idea what you intended to do or why the code is now failing, other than I'm pretty sure zed[[[2]]] is part of the cause. Actually I think you might have edited and changed your question since since I posted my proposed answer. $\endgroup$ – Bill Jul 11 '17 at 23:58

Your Answer

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

Browse other questions tagged or ask your own question.