4
$\begingroup$

I have a set of three coupled differential equations and want to explore the properties of the system by applying various techniques in the field of dynamical systems. I found a package here, called EcoEvo, which seems it has built-in functions for various analyses in this regard; but I cannot make it work. First, one needs to install the package:

PacletInstall["EcoEvo", "Site" -> "http://raw.githubusercontent.com/cklausme/EcoEvo/master"]

Then, to load it and to set up my model, I write:

<< EcoEvo`

SetModel[{Pop[p] -> {Equation :> -gamma p[t] f[t]}, Pop[f] -> {Equation :> -c f[t] + gamma p[t] f[t]}, Pop[T] -> {Equation :> b (1 - T[t]) (1 - 0.214499/T[t]) - m (1/(1 + (r/f[t])^n)) (1 - 0.214499/T[t])}}]

At this stage, an error is returned, but it seems it is not problematic and my system is defined. Now, for example, let's find the equilibria for a given set of parameters:

gamma = 0.552; c = 0.004; b = 0.093; r = 0.798; n = 10.862; m = 0.870;
eq = SolveEcoEq[]

Nothing is returned. Am I executing the commands correctly?

$\endgroup$
7
  • 3
    $\begingroup$ It's likely that you will have to contact the package developer for better support. This does not seem a widespread package, although diffeq is not my field and I may simply not be familiar with it. Having said that, what error do you receive at the start? It may help to mention that explicitly. The GitHub page mentions that documentation is available both online and in the package, but unfortunately the online one is not accessible. Have you tried to run a sample problem from the documentation to check that your setup is correct? $\endgroup$
    – MarcoB
    Commented Feb 14 at 11:17
  • 3
    $\begingroup$ Note that according to the documentation, the variables do not have the time included. Namely, don't write p[t], f[t] and T[t], but just p and f and T. Second, it is not that nothing is returned, what is returned is an empty set. This indicates that there is no equilibrium. You can check that also manually with Solve: sol = Solve[{0 == -gamma p f, 0 == -c f + gamma p f, 0 == b (1 - T) (1 - k/T) - m (1/(1 + (r/f)^n)) (1 - k/T)}, {p, f, T}, MaxExtraConditions -> All] None of the solutions returned is compatible with your parameters. $\endgroup$
    – Domen
    Commented Feb 14 at 11:38
  • 2
    $\begingroup$ @Anovice, that's really not the case. Think about a simple ODE such as $y'=1$. Obviously, it doesn't have any equilibrium. $\endgroup$
    – Domen
    Commented Feb 14 at 13:12
  • 1
    $\begingroup$ I've made a small mistake; $p$ is not necessarily zero, it can be any number, so $(p,0,1)$ and $(p,0,k)$ (but the issue with $f\to 0$ remains). $\endgroup$
    – Domen
    Commented Feb 14 at 14:56
  • 3
    $\begingroup$ @MarcoB Oops, the online documentation should be back now. It's currently hosted on Wolfram Cloud, which helpfully deletes my files every 60 days. I hope to move to the Wolfram Paclet Repository in the near future. $\endgroup$
    – Chris K
    Commented Feb 15 at 0:55

1 Answer 1

4
$\begingroup$

As @Domen noted, the package doesn't use [t] in defining a model since v1.6.0. At some point I'll fix old answers on this site to reflect that change.

On to your question. Multiplying the last term of T'[t] by f^n/f^n gets rid of one possible divide-by-zero situation, which helps get some results.

<< EcoEvo`

SetModel[{
  Pop[p] -> {Equation :> -gamma  p  f},
  Pop[f] -> {Equation :> -c  f + gamma  p  f}, 
  Pop[T] -> {Equation :> b  (1 - T)  (1 - 0.214499/T) - m (f^n/(f^n + r^n))  (1 - 0.214499/T)},
  Parameters :> {gamma, c, b, m, r, n}
}]

gamma = 0.552; c = 0.004; b = 0.093; r = 0.798; n = 10.862; m = 0.870;

eq = SolveEcoEq[]

(* some warnings, including *)
(* Solve::svars Equations may not give solutions for all "solve" variables. *)
(* {{f -> 0, T -> 0.214499}, {f -> 0, T -> 1.}} *)

This is because if $f=0$, then $p$ can be anything -- a bit degenerate. But $p=0$ seems like the best bet since $p$ approaches zero. The following plots a 2D phase-plane between $f$ and $T$ assuming $p=0$.

PlotEcoPhasePlane[{f, 0, 1}, {T, 10^-10, 2}, Fixed -> {p -> 0}]

enter image description here

This shows the two equilibria. A numerical solution verifies that $f=0, T=1$ is stable.

sol = EcoSim[{p -> 1, f -> 1, T -> 2}, 1000];
PlotDynamics[sol]

enter image description here

You can see how $T$ follows its isocline as $f\to 0$.

If you start below the $T$ isocline, the system diverges and you run into another divide-by-zero situation when $T$ hits zero.

sol = EcoSim[{p -> 1, f -> 1, T -> 0.1}, 1000];
PlotDynamics[sol]
(* NDSolve::ndsz At t == 370.58301686040124`, step size is effectively zero; singularity or stiff system suspected. *)

enter image description here

If $T$ has some physical interpretation as a non-negative quantity, you might rethink your equations, because they don't respect that.


Using OP's initial conditions illustrates that you need to be careful with AccuracyGoal. Using the default:

sol = EcoSim[{p -> 0.817, f -> 0.001, T -> 0.85}, 500];
PlotDynamics[sol]

enter image description here

🪦 RIP

Using more accurate solution:

sol = EcoSim[{p -> 0.817, f -> 0.001, T -> 0.85}, 500, NDSolveOpts -> {AccuracyGoal -> \[Infinity]}];
PlotDynamics[sol]

enter image description here

❤️‍🩹 Recovery!

$\endgroup$
5
  • $\begingroup$ While the $(1−0.214499/T)$ term vanishes at $T=0.214499$, if $T<0.214499$ it's negative and if $T=0$, $dT/dt$ becomes infinitely negative which causes numerical problems. If the physical meaning of $T$ is such that it should remain positive, this is also a conceptual problem with your model. What do the variables represent? $\endgroup$
    – Chris K
    Commented Feb 15 at 18:46
  • $\begingroup$ I know nothing about the system, but here are my naive impressions. I'm not sure the difference between $f$ and $T$ -- should $f$ actually be a variable or just have $p$ affect $T$? I'd also want more justification of the functional form of interactions (e.g. why does $f$ affect the rate $p$ changes, why $1/T$, etc.) Where are the parameters from? This chapter from Ellner & Guckenheimer might provide some modeling advice. Good luck! $\endgroup$
    – Chris K
    Commented Feb 16 at 0:42
  • $\begingroup$ @Anovice The phase-plane is 2D (assuming $p=0$) so it's not clear how to project the 3D system onto it. Anyhow, $T$ can become negative in the phase-plane because it can in the differential equations for the right set of ICs. $\endgroup$
    – Chris K
    Commented Feb 16 at 17:59
  • $\begingroup$ @Anovice To me, the point of a dynamic model is to include the mechanisms at work. Then the model can not only recreate the phenomenon you saw but also extrapolate to other conditions. So if it produces physically unrealistic results for some realistic initial conditions, I'd be concerned. But your needs may differ. $\endgroup$
    – Chris K
    Commented Feb 17 at 16:44
  • 1
    $\begingroup$ No, what happens in one part of phase space is not necessarily related to what happens in another part. I can't exactly say your model has a stable equilibrium, because it is degenerate in another way: there is a line of equilibria where $f=0, T=1$, and $p$ can be anything. The eigenvalues of the Jacobian matrix at $f=0, T=1, p=0$ for example are {-0.0730516, -0.004, 0}, which is consistent with this neutral stability. If you want to get more background on nonlinear differential equations, try the book "Nonlinear dynamics and Chaos" by Steven Strogatz. $\endgroup$
    – Chris K
    Commented Feb 27 at 2:59

Your Answer

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