3
$\begingroup$

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?

Improve this question
$\endgroup$
4
  • $\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
    Commented Oct 7, 2014 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$
    – Ramiro Magno
    Commented Oct 7, 2014 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
    Commented Oct 8, 2014 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
    Commented Oct 8, 2014 at 2:59

1 Answer 1

Reset to default
1
$\begingroup$

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]
]}}
Improve this answer
$\endgroup$
2
  • $\begingroup$ Indeed, would you know how I could use that solution for fitting instead? $\endgroup$
    – Ramiro Magno
    Commented Oct 17, 2014 at 5:48
  • $\begingroup$ Sorry I don't have much experience with fitting models. Maybe someone else could chime in? $\endgroup$
    – Chris K
    Commented Oct 18, 2014 at 15:38

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.