I am trying to use this model for later use in FindFit:

model[x0_?NumberQ, r_?NumberQ, k_?NumberQ, n_?NumberQ] := model[x0, r, k, n] = 
   First[x /. NDSolve[{x'[t] == r x[t] (1 - x[t]/k)^n, x[0] == x0}, x, {t, 0, 61}]]

Assume now that those parameters are sought within the following ranges:

x0: 10^-50..10^-10;
n: 0.1..5.0;
r: 1..10;
k: 1..1000.

Depending on the specific choice of the parameter combination I will get various NDSolve numeric errors. So what would be a good choice of MachinePrecision, AccuracyGoal and PrecisionGoal that would work for that parameter space?

  • $\begingroup$ I have an idea what the issue is, but I wanted to be sure I understand what you're asking. I've tried several inputs and can get no error messages (V10.0.1, Mac OSX). Or do you mean only that the computed result looks wrong? $\endgroup$ – Michael E2 Oct 7 '14 at 13:34
  • $\begingroup$ It varies. Can be NDSolve error messages such as NDSolve::precw: The precision of the differential equation (...) is less than WorkingPrecision (40.), or indeed simply the InterpolatingFunction returned (without NDSolve errors) is just nonsensical. $\endgroup$ – plant Oct 7 '14 at 14:06
  • $\begingroup$ I won't have time to write up an answer for a few days, but here are couple things to consider. For n >= 1, the options PrecisionGoal -> 10, AccuracyGoal -> Infinity, WorkingPrecision -> 100 seem to take care of everything I tried. A WorkingPrecision around 100 seems to be needed for x0 as small as 10^-50 -- the smaller x0, the greater the precision needed. Second, I think for n < 1, x reaches k in finite time. I haven't figured out a graceful way to stop the integration at the right time. Also n = 1.01 behaves like n < 1 -- not sure why exactly. $\endgroup$ – Michael E2 Oct 8 '14 at 2:57
  • $\begingroup$ The message NDSolve::precw can be ignored. Ideally you would pass parameters x0, r, k, and n with a precision at least as great as the WorkingPrecision -- then you should not get the message. $\endgroup$ – Michael E2 Oct 8 '14 at 2:59

This looks like a variant of the logistic growth population model, so I'll call x "population size". A few thoughts:

First, in such cases I've had good luck with AccuracyGoal->Infinity when dealing with low population sizes. This is because a small absolute error in population size will be critical when you're around the unstable equilibrium x=0. This seemed to fix any problem I had with your model (as long as n wasn't too small).

Second, it might also help to deal with log population sizes, since exponential population growth is linear on a log scale. That is,

model2[x0_?NumberQ, r_?NumberQ, k_?NumberQ, n_?NumberQ] := model2[x0, r, k, n] =
First[lnx /. NDSolve[{lnx'[t] == r  (1 - E^lnx[t]/k)^n, lnx[0] == Log[x0]}, lnx, {t, 0, 61}]]

Third, it seems like this particular model has an analytical solution available through DSolve:

DSolve[{x'[t] == r x[t] (1 - x[t]/k)^n, x[0] == x0}, x, t]


{{x -> Function[{t}, InverseFunction[
-((Hypergeometric2F1[n,n,1+n,k/#1](1-k/#1)^n (1-#1/k)^-n)/n)&][rt-((1-k/x0)^n (1-x0/k)^-n Hypergeometric2F1[n,n,1+n,k/x0])/n]
  • $\begingroup$ Indeed, would you know how I could use that solution for fitting instead? $\endgroup$ – plant Oct 17 '14 at 5:48
  • $\begingroup$ Sorry I don't have much experience with fitting models. Maybe someone else could chime in? $\endgroup$ – Chris K Oct 18 '14 at 15:38

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.