# 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 == 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? – Michael E2 Oct 7 '14 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. – rmagno Oct 7 '14 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. – Michael E2 Oct 8 '14 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. – 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 == 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 == 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? – rmagno Oct 17 '14 at 5:48
• Sorry I don't have much experience with fitting models. Maybe someone else could chime in? – Chris K Oct 18 '14 at 15:38