3

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?

Sektor
  • 3,320
  • 7
  • 27
  • 36
Ramiro Magno
  • 325
  • 1
  • 10
  • 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 07 '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. – Ramiro Magno Oct 07 '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 08 '14 at 02: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 08 '14 at 02:59

1 Answers1

1

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]
]}}
Chris K
  • 20,207
  • 3
  • 39
  • 74