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.
 A: 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.


A: 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$.
