Transition time in a Lotka-Volterra system

Question

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.