2
$\begingroup$

I tried to search around this website for an answer and found many similar requests, but for some reason I cannot get any of them to work.

Anyway, I am trying to fit data to a solution of 3 ODE's. The data is something like this.

data = {{0, 0}, {1, 0}, {2, 0}, {3, 139}, {4, 261}, {5, 452}, {6, 623}, {7, 
806}, {8, 992}, {9, 1189}, etc etc }

The code I am using is the following. Just to note, these ODE's plot just fine when I use NDSolve and set k1,k2,k3 to particular values. So I am assuming there is something wrong which I am ignorant to. I am not super familiar with this. I tried to follow this to no avail.

Also, it just keeps working. I do not seem to get any errors but I even let it run overnight to no avail. (Sorry for using a jpg, I failed to get the code to look correct in the editor). The full data set is here, https://github.com/tpkdesigns/deleted.git if anyone needs it.

enter image description here

I appreciate anyone's time that is put into this.

Edit Per request here is the actual code. It copy and pastes strangely(greek letter and derivatives), that's why I didn't do it originally.

Clear["Global`*"]
currentdir = "/Users/Physics/Desktop/";
SetDirectory[currentdir];
\[Gamma] = 6;
data0det = Import["deleted.txt", "Table", HeaderLines -> 1];
data = data0det[[All, {1, 4}]];
ode = { Derivative[1][Nb][t] == 
    k1*Na[t] - \[Gamma] Nb[t] + k2*Na[t]*Nc[t]*Nc[t],
   Derivative[1][Na][t] == -k1* Na[t] + 1/2 \[Gamma] Nb[t],
   Derivative[1][Nc][t] == 1/2 \[Gamma] Nb[t] - k2*Na[t]*Nc[t]*Nc[t],
   Na[0] == k3, Nb[0] == 0, Nc[0] == 0 };
model = ParametricNDSolveValue[ode, Nb, {t, 0, 956}, {k1, k2, k3}]
fit = NonlinearModelFit[data, 
  model[k1, k2, k3][t], {{k1, 10}, {k2, 20}, {k3, 30}}, t]
fit["BestFitParameters"]
plotfit = Plot[fit[t], {t, 0, 200}]
$\endgroup$
  • 1
    $\begingroup$ Please do not post images of your work, especially when the images display at a size that make them difficult to read. Please post your actual Mathematica code in the form of text that can be copied and pasted into a Mathematica notebook. Without such, it will be difficult to reproduce your problem and to experiment with possible solutions. $\endgroup$ – m_goldberg May 4 '18 at 13:29
  • $\begingroup$ @m_goldberg Code posted sir. $\endgroup$ – B.Jones May 4 '18 at 13:41
  • 1
    $\begingroup$ Should $k_1,k_2,k_3$ be positive? or are these parameters allowed to be negative? $\endgroup$ – AccidentalFourierTransform May 4 '18 at 14:36
  • $\begingroup$ @AccidentalFourierTransform they should be positive. k1 is a pump rate, k2 is a depump rate, and k3 is just some positive number. $\endgroup$ – B.Jones May 4 '18 at 14:39
  • $\begingroup$ @B.Jones Your starting parameters do not represent the shape of your data at all. Your data contains a peak, whereas your model with those values of $k_i$ just rises to a plateau. Can you find better starting points? The scale of the model's response vs. the data also seems way off, but a few orders of magnitude. $\endgroup$ – MarcoB May 4 '18 at 15:50
1
$\begingroup$

With $\gamma=6$ it seems impossible to fit your data accurately. So I'm going to include this parameter as one of the free variables. Also, it seems that the best fit requires $k_2=0$, so I'm going to set this parameter to zero from the beginning to accelerate the convergence. With this,

model[γ_?NumericQ, k1_?NumericQ, k2_?NumericQ, k3_?NumericQ] :=
model[γ, k1, k2, k3] = NDSolveValue[{nb'[t] == 
  k1*na[t] - γ nb[t] + k2*na[t]*nc[t]*nc[t], 
  na'[t] == -k1*na[t] + 1/2 γ nb[t], 
  nc'[t] == 1/2 γ nb[t] - k2*na[t]*nc[t]*nc[t], na[0] == k3,
  nb[0] == 0, nc[0] == 0}, {na, nb, nc}, {t, 0, 956}][[2]]

NonlinearModelFit[data, model[γ, k1, 0, k3][t], {{γ, .0399887}, {k1, .0560411}, {k3, 4000.5}}, t]

Show[Plot[%[t], {t, 0, 950}, PlotStyle -> Red, PlotRange -> All], ListPlot[data]]

enter image description here

The fit is not perfect, but I'm not sure whether it's because of a poorly chosen starting point, or because the model itself is not a good representation of the system. I leave to you to explore this.

$\endgroup$
  • $\begingroup$ FWIW: the data seems to follow a sort of generalised Cauchy distribution, $\frac{a}{b+x^{2/5}}$, or an inverse gamma distribution, $x^{-a}\exp(-b/x)$ (with suitable $a,b$, and shifting, if necessary, the mean $x\to x-\mu$). $\endgroup$ – AccidentalFourierTransform May 4 '18 at 22:53
  • $\begingroup$ I really appreciate you for looking into this. Yes it does seem to work now pretty well but you are right that the model may not represent the physics of the system. That conclusion is fine with me and I will have to explore the details further. $\endgroup$ – B.Jones May 8 '18 at 13:14
  • $\begingroup$ @B.Jones I'm glad I could help :-) If you need any more help, feel free to ask. Cheers! $\endgroup$ – AccidentalFourierTransform May 8 '18 at 13:19

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.