Asymptotic solution of a system of ODEs

I have asked this question on math.stackexchange, however, have not got any answer. Therefore, I suspect that this system of ordinary differential equations cannot be solved analytically. But I still have some hope that it is possible to show that $$a(u)\sim -1/u$$ for $$u\rightarrow \infty$$. The system arises in physics, approximate treatment of the Hubbard model using some newly developed formalism.

$$2 a'(u)=u x'(u),$$ $$4 a u x'(u)-4 y'(u)=-2 x(u) \big[u a'(u)+a(u)\big],$$ $$4 a u y'(u)-4 x'(u)=\big[u\lambda'(a)a'(u)+\lambda(a)\big]-2y(u)\big[u a'(u)+a(u)\big].$$ with initial conditions $$x(0)=y(0)=0, a(0)=-1/2,$$ and $$\lambda(a)=a^3+\frac14 a$$.

It is easy to find a series solution in the vicinity of $$u=0$$. It is also possible to integrate one of the equations. But the problem is in the infinity point $$u\rightarrow \infty$$.

To my astonishment, a simpler approximation results in a system of 2 ODEs of similar structure, which indeed permits analytic solution.

Edit

In the meantime I managed to integrate 2 equations out of 3 and got: $$y(u)=\frac{1}{2} u a(u) x(u)+\frac{a(u)^2}{2}-\frac{1}{8},\\ 2 \Lambda (a(u))-u x(u) \lambda (a(u))+4 u a(u) x(u) y(u)-2 \left(x(u)^2+y(u)^2\right)=\frac{3}{32},$$ where $$\Lambda (a)=\frac{1}{8} a^2 \left(2 a^2+1\right)$$. But I still do not know how to get asymptotics. The remaining equations are $$2 a'(u)=u x'(u),\\ 4 a(u)^2 \left(3 u^2 x(u)^2+1\right)=4 x(u) \big(u a(u)+4 x(u)\big)+1.$$

• Numerically, the large-$u$ behavior is given by $a(u)= -1.02418609585 u^{-1} +0.56465348534 u^{-3} +1.36082946003 u^{-5} +0.94254923789 u^{-7}$, $x(u)= 0.34089869589 -1.02418609585 u^{-2} +0.84698022802 u^{-4} +2.26804910005 u^{-6}$, $y(u)= -0.2995718522 +1.1452019773 u^{-2} -1.0692455654 u^{-4} -2.6928672430 u^{-6}$ for your initial conditions. (Mathematica, high-precision NDSolve, expansion around $u=\infty$) Commented Jan 10 at 14:24
• @FredHucht Yes, this is quite plausible. Commented Jan 10 at 14:27
• @FredHucht Why not posting your solution? At first I was convinced that $a(u)=-1/u +\ldots$ and that deviation in the leading coefficient of your solution is just a numerical noise. But I am more convinced now that it is real. Commented Jan 12 at 19:06
• I'm quite sure that the numerics are correct. The leading term -1.024... varies continuously with $a(0)$. I do not find a better way to handle your ODE system. One could use your edit and write down one single ODE, but this is highly nonlinear. Maybe I play around with it today, then I'll write an answer. Commented Jan 12 at 19:24
• Commented Jan 21 at 7:12

This is my comment as an answer. Using Mathematica's ParametricNDSolve[] with 50 digits precision as well as a series expansion around $$u=\infty$$, I get \begin{align} a(u) &= -1.02418609585 u^{-1} + 0.56465348534 u^{-3} + O(u^{-5}) \\ x(u) &= 0.34089869589 - 1.02418609585 u^{-2} + O(u^{-4}) \\ y(u) &= -0.2995718522 + 1.1452019773 u^{-2} + O(u^{-4}) \end{align} for the given initial condition $$a_0=a(0)=-1/2$$. Therefore, the leading term $$a_\infty=\lim_{u\to\infty} u\,a(u)$$ is $$a_\infty=-1.024\ldots \neq -1$$. We can however ask the inverse question, namely for which $$a_0^*$$ we get $$a_\infty=-1$$, with result $$a_0^*=-0.92708111527$$.
Interestingly, the Taylor series of $$a(u)$$ around $$u=\infty$$ contains a lot of terms proportional to $$\sqrt{a_\infty^2-1}$$, and significantly simplifies for $$a_\infty=-1$$. As example, the $$u^{-3}$$ term reads $$\frac{a_{\infty } (4-3 a_{\infty }^2) }{3 a_{\infty }^2+1} \left(3 a_{\infty } \sqrt{a_{\infty }^2-1}-2\right)\stackrel{a_\infty=-1}{\mapsto} \frac 1 2 \,.$$ Therefore, the limit $$a_\infty=-1$$ seems to be a special case.
Similarly, the expansion around $$u=0$$ contains terms proportional to $$\sqrt{4a_0^2-1}$$, which drop out for $$a_0=-1/2$$.
• Thank you, very useful! In other words, the $u^{-1}$ coefficient can only be obtained numerically and is different from $-1$. On the other hand, all other coefficients of the series expansion can be written in terms of it. Commented Jan 13 at 22:59
• Yes, and it depends on $a_0$ (assuming $x(0)=y(0)=0$). Maybe, $a_\infty=-1$ holds if you slightly change the ODE and especially $\lambda(a)$. Commented Jan 13 at 23:05