# Reasonable choice of MachinePrecision, AccuracyGoal and PrecisionGoal for NDSolve

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?

• 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? Commented Oct 7, 2014 at 13:34
• 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. Commented Oct 7, 2014 at 14:06
• 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. Commented Oct 8, 2014 at 2:57
• 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. Commented Oct 8, 2014 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]


gives

{{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]
]}}

• Indeed, would you know how I could use that solution for fitting instead? Commented Oct 17, 2014 at 5:48
• Sorry I don't have much experience with fitting models. Maybe someone else could chime in? Commented Oct 18, 2014 at 15:38