Transition time in a Lotka-Volterra system

I am working with a set of real-valued ordinary differential equations based on the Lotka-Volterra competition equations:

\begin{align} \dot{a_1} & = a_1 \left( 1 - a_1 - 2 a_2 \right) \\ \dot{a_2} & = a_2 \left( 1 - a_2 - (1-1/\nu) a_1 \right) \end{align}

where $a_{1,2} \in [0,1]$ and $\nu \ge 1$. I would like to obtain a closed form (or analytical) solution for the time, $\tau$, it takes for this system to transit between two regions in state space. Specifically, I would like to solve for $\tau$ given $a_1(0) = \delta$ and $a_2(\tau) = \delta$, where $1/2 \lt \delta \lt 1$. This is along the manifold between two fixed points: the unstable manifold of a saddle at $(a_1,a_2) = (1,0)$ and a stable equilibrium at $(a_1,a_2) = (0,1)$.

There is a solution for $\nu = 1$ ($\tau = -2 \text{ln}(1/\delta-1)$), but I have been unable to find a solution for the more general case in terms of $\nu$ (or even for any particular value of $\nu \gt 1$, including the limiting case of $\nu \rightarrow \infty$). This paper seems to imply that the equations above are not fully integrable except when $\nu = 1$ (see Eq. 25'), i.e, when the second equation is not a function of $a_1$. However, I'm not actually interested in solving this system for all time over the full state space. Question: Are there any methods to solve or obtain a reliable approximation for $\tau = f(\nu, \delta)$ for this system just within my region of interest?

Attempts:

• In addition to paper and pencil, I have used Matlab's dsolve and Mathematica's DSolve along with assumptions to try to solve the ODEs for the specified boundary conditions. I was unable solve the system using these, but might there be ways to transform or break up the system that would facilitate a solution?

• I have tried using low-order power series, e.g., this paper, about each of the equilibria to obtain functions of time that are then inverted to solve for $\tau$. This was far from accurate as only a few series terms can be used.

• I have tried schemes based on simulating the ODE and fitting the transit times to a function of $\nu$. This requires finding initial conditions that lie on the manifold. How can that be done reliably as a function of $\nu$ and $\delta$? And can fitting methods be adapted to more general forms of the equation above (i.e., reusing the same fitting coefficients and without fitting high dimensional surfaces)? I am more interested in non-fitting-based solutions, but if you can demonstrate something that works well, I would be happy to look at it.

Update 1 – Mar. 17, 2014:

The nullclines of the system and the Jacobian determinant (related to curvature) can be used to obtain estimated initial and final conditions. Solving for the roots of $\det(J)=0$ as a function of $a_1$ and scaling appropriately, one obtains

$$a_2(t) = \tfrac{1}{2} \left( 2 - a_1(t) - \sqrt{a_1(t) \left( \tfrac{2}{\nu} \left( a_1(t)-1 \right) + 2 - a_1(t) \right)} \right)$$

For $a_1(0) = 1-\delta$, this expression appears to provide a good approximation for $a_2(0)$ on the stable manifold (attracting contour) in question. It is less reliable for obtaining an estimate for $a_1(\tau)$. Perhaps this is due to the Jacobian evaluated at $(a_1,a_2) = (0,1)$ being a defective matrix with only one eigenvector.

Update 2 – Mar. 25, 2014 – (post bounty):

The equation from @JJacquelin's answer immediately after the transformation to polar coordinates can be simplified and integrated definitely (I used Mathematica 9) with respect to the bounds at times $0$ and $\tau$:

\begin{align} \int_0^\tau{dt} & = \int_{\theta_0}^{\theta_{\tau}}{\frac{\sin^3\theta + \cos^3\theta + \left(2\cos\theta + \left(1-\tfrac{1}{\nu}\right)\sin\theta\right)\sin\theta\cos\theta}{\left(\tfrac{1}{\nu} \cos \theta + \sin \theta\right)\sin\theta\cos \theta}d\theta} + \int_{\rho_0}^{\rho_{\tau}}{\frac{1}{\rho}d\rho} \\ \tau & = \int_{\theta_0}^{\theta_{\tau}}{\frac{1-\tfrac{1}{\nu}+\left(2+\cot\theta\right)\cot\theta+\tan\theta}{1+\tfrac{1}{\nu}\cot\theta}d\theta} + \ln\left(\frac{\rho_{\tau}}{\rho_0}\right) \\ & = \frac{\nu^2}{1+\nu^2}\left(\ln\left(\frac{\cos(\theta_0)}{\cos(\theta_{\tau})}\right) +2\ln\left(\frac{\cos\theta_{\tau}+\nu\sin\theta_{\tau}}{cos\theta_0+\nu\sin\theta_0}\right)+\frac{1}{\nu^2}\ln\left(\frac{\cos\theta_{\tau}\left(1+\nu\tan\theta_{\tau}\right)^2}{\cos\theta_0\left(1+\nu\tan\theta_0\right)^2}\right)\right) + \nu\ln\left(\frac{\nu+\cot\theta_0}{\nu+\cot\theta_{\tau}}\right) + \ln\left(\frac{\rho_{\tau}}{\rho_0}\right) \end{align}

Transforming back from polar to Cartesian coordinates further simplifies the expression:

$$\tau = \ln\left(\frac{a_1(0)}{a_1(\tau)}\right) + 2\ln\left(\frac{a_1(\tau)+\nu a_2(\tau)}{a_1(0)+\nu a_2(0)}\right) + \nu\ln\left(\frac{a_2(\tau)\left(a_1(0)+\nu a_2(0)\right)}{a_2(0)\left(a_1(\tau)+\nu a_2(\tau)\right)}\right)$$

Here is some "simple" Matlab code that demonstrates this. Also in the code is an ODE-based method that uses ode45's event detection. This accurately finds the time $\tau$ by integrating along the manifold in question and terminating when the solution satisfies a condition. It's simple and fast for this basic case. However, recall that this system is a simplified version of a more general one in which the time scaling of the vector field can vary and even be different in each of the two dimensions. This could lead to long integration times if the scaling cannot be factored. I'm hoping that the analytical solution given here can be generalized.

The analytic solution-based part of my Matlab code underestimates $\tau$ (except for $\nu$ close to $1$ when it slightly overestimates it). The source of almost all of the error is due to the estimate of $a_1(\tau)$ (af1 in the code). Something else I'm looking into.

• "For a1(0)=1−δ, this expression appears to provide a good approximation for a2(0) on the stable manifold" Could you expand on this? Are these points supposed to be on the "manifold"? – Did Mar 17 '14 at 7:03
• @Did: I can't prove that the contour of zero curvature corresponds to the manifold (contour) along which I trying to find the time. It's very very close, but but numerical simulations seem to show that there's a difference, especially as $a_1(0)=\delta\rightarrow1/2$. I haven't had the chance to investigate this closely. – horchler Mar 17 '14 at 17:51
• could you please send me the Matlab codes ? Thanks – user168724 Aug 8 '14 at 19:41
• @baouche: Please be specific about which Matlab code you mean. I provided a link to my own code in the second update to my question. – horchler Aug 8 '14 at 21:09

Knowing the analytical solution might be of interest for your work, even if the solution is only provided on a parametric form in polar coordinates.

The general solution above is given without initial conditions. If the initial conditions $a_1(0)$ and $a_2(0)$ are given, the problem is fully determined and a unique solution can be derived :

Outside of the question initially strictly raised, a general view of the solutions of the system of two ODEs is shown below : The map of $a_1 , a_2$ and a sketch of the behaviour.

• Interesting. However, this is in terms of indefinite integrals. It's not clear to me that explicit expressions for the constants of integration can be obtained as a function of $\delta$ (or even some $\theta_0$). I need too spend some time with this, but clarify this that would be helpful. – horchler Mar 14 '14 at 20:17
• In order to explicit the constants of integration, the initial conditions must be defined (as updated above). Nevertheles, in order to study the functions $a_1(t) , a_2(t)$ and to compute the value of $t$ to go from the initial point to the final point, the use of the analytic formulas is not the best method : the formulas are complicated and it is not an easy task to write a program without mistake. It is much easier to solve numerically the system of two EDOs. This requires only a few lines of code which can be written in a few minutes. – JJacquelin Mar 17 '14 at 9:49
• A few seconds in Matlab actually. I know what you mean. This is more or less what I do now, but I'd rather not perform numerical simulation of the ODEs for every instance of $\nu$ and $\delta$ (and other parameter in the more general case). This is part of an optimization problem so the integration might have to be performed on every iteration. By the way, have you noticed that your $f(\theta)$ integral can be "simplified" in terms of the Gaussian hypergeometric function $_2F_1$? Usually an "analytic solution" wouldn't contain integrals. – horchler Mar 17 '14 at 16:13
• I agree with your approach. Of course, Matlab can do the simulations as well (and probably better) than a simplified program such as I suggested. That is great that the integral can be simplified on the form of a Gauss hypergeometric function. I even not seach for a closed form of the integral because I think that it is not useful in practice, as I already said. But on a pure theoretical viewpoint, it is of interest. – JJacquelin Mar 17 '14 at 17:08
• Your integrand of $f(\theta)$ at the end seems to be $1/g(\theta)$, which doesn't seem to match your earlier equations. Typo somewhere? – horchler Mar 17 '14 at 17:09

We first define rigorously the set $S$ of interest, named "the stable manifold" in the question. To do so, note that the nonnegative quadrant $x\geqslant0$, $y\geqslant0$, is stable by the dynamics, and that, for every starting point in the positive quadrant $x\gt0$, $y\gt0$, the dynamics converges to $(0,1)$. Consider the reversed dynamics, defined as $$x'=-x\cdot(1-x-2y),\qquad y'=-y\cdot(1-y-(1-1/\nu)x).$$ Every solution of the reversed dynamics starting from $(x_0,y_0)$ in the strip $x\in(0,1)$, $y\gt0$, either reaches the line $x=1$ at some point $(1,y)$ with $y\gt0$, or at $(1,0)$, or converges to $(0,0)$.

The "stable manifold" from $(1,0)$ to $(0,1)$ for the original dynamics is the set $S$ of initial conditions $x_0\in(0,1)$, $y_0\gt0$, such that the solution of the reversed dynamics reaches $x=1$ at $(1,0)$. Very loosely speaking, $S$ can be viewed as the trajectory of the original dynamics starting from $(1,0)$.

Second we revisit the case $\nu=1$. Then, $S$ is the line $a_1+a_2=1$ hence $\tau$ is the time to go from $(a_1,a_2)=(\delta,1-\delta)$ to $(a_1,a_2)=(1-\delta,\delta)$. The differential equation $a_2'=a_2(1-a_2)$ has an explicit solution, namely, $$\frac{a_2(t)}{1-a_2(t)}=\mathrm e^t\cdot\frac{a_2(0)}{1-a_2(0)}.$$ Plugging in this identity the values $a_2(0)=1-\delta$ and $a_2(\tau)=\delta$ yields $$\tau=2\log\left(\frac{\delta}{1-\delta}\right),$$ that is, the formula in the post (with a minus sign missing). More generally, the time to go from $a_2(0)$ to $a_2(\tau)\geqslant a_2(0)$ is $$\tau=\log\left(\frac{a_2(\tau)\cdot(1-a_2(0))}{(1-a_2(\tau))\cdot a_2(0)}\right),$$ hence, when $a_2(0)\to0$ and $a_2(\tau)\to1$, $$\tau\sim-\log a_2(0)-\log(1-a_2(\tau)).$$

Consider finally $\nu\gt1$. Let $\mu=1/\nu$. On the line $(1-\mu)a_1+a_2=1$, $a_2'=0$ and $a'_1\lt0$. On the line $a_1+2a_2=1$, $a_1'=0$ and $a'_2\gt0$. Thus the region between these two lines is stable and $S$ lies entirely between these two lines. This implies that, everywhere on $S$, $$\mu a_2(1-a_2)\leqslant a_2'\leqslant a_2(1-a_2).$$ Considering the differential system near the point $(a_1,a_2)=(1,0)$, one sees that, when $\delta\to1$, the point on $S$ such that $a_1=\delta$ is such that $a_2\sim(1-\delta)^\mu$, thus the time $\tau$ to go from $a_1=\delta$ to $a_2=\delta$ is roughly the time to go from $a_2=(1-\delta)^\mu$ to $a_2=\delta$. The exact result in the case $\nu=1$ and the double inequality on $a_2'$ yield $$\liminf_{\delta\to1}\left(\frac{\tau}{-\log(1-\delta)}\right)\geqslant\mu(1+\mu),$$ and $$\limsup_{\delta\to1}\left(\frac{\tau}{-\log(1-\delta)}\right)\leqslant1+\mu.$$ Loosely speaking, for some $\mu(1+\mu)\leqslant\alpha\leqslant1+\mu$, $$\mathrm e^\tau\approx(1-\delta)^{-\alpha}.$$ Since $\alpha\lt2$, this is a speed-up when compared to the $\nu=1$ case, due to the fact that starting from $a_1(0)=\delta$ on $S$ is to start from $a_2(0)\gg1-\delta$.

• I like the idea of bounding the time. It would be useful to obtain $\tau$ for the limiting case of $\nu \rightarrow \infty$. Some questions. Are the two bounding lines chosen arbitrarily – could tighter bounds be used? Why not $a_1+a_2 = 1$ for the upper bound? – horchler Mar 16 '14 at 17:51
• "Why not $a_1+a_2=1$ for the upper bound?" Hmmm... the line $a_1+a_2=1$ is below the manifold, at least for $a_1\sim1$ (and probably in its entirety). – Did Mar 16 '14 at 18:04
• We may be using "manifold" for different things. I'm not a mathematician. I'm referring to the unique contour that all trajectories are attracted to. I'm only interested in the time along this contour (or finding initial conditions on it). – horchler Mar 16 '14 at 18:10
• Trajectories are not attracted to anything except to the point $(0,1)$... but anyway, I am following you in the meaning of "manifold". A neat definition of this set would be as the set of starting points such that the trajectory of the reversed dynamical system hits the line $a_1=1$ at $a_2=0$ (above it one hits the right side of the square at some point $(1,a_2)$ with $a_2\gt0$, and below it, one goes to $(0,0)$). – Did Mar 16 '14 at 18:16
• I've been trying to extract useful bounds from this. Your estimate of $\tau\approx(1-\delta)^{-\alpha}$ doesn't seem to work for small $\nu$. In particular it greatly overestimates the case of $\nu=1$. In simulations, I find that the time to go from $a_2(0)=(1-\delta)^\mu$ to $a_2(\tau)=\delta$ underestimates the time from $a_1(0)=\delta$ to $a_2(\tau)=\delta$ (moreso for large $\nu$). I'm not sure how I can make practical use of these results. – horchler Mar 16 '14 at 20:56