Analytic solution expected

Question

I am desperately trying to solve a system of differential equations (post at mathoverflow helped me to clarify the asymptotic properties of the solution). In the meantime I've reduced 3 ODEs to just one

Clear[p,q,r]
p[a_,u_]:=-(16+u^2)^2+3 u^2 a^2 (36 u^4*a^4+3*u^2*(u^2-48)*a^2-u^4+8*u^2+192)
q[a_,u_]:= 2*u*a*(9*u^4*a^4-3*u^2*(u^2-4)*a^2-(u^2+16))
r[a_,u_]:= u^2*a^2*(3*u^2*a^2*(3*a^2-1)-1)
eq=r[a[u],u]+q[a[u],u]*a'[u]+p[a[u],u]*a'[u]^2==0
ic=a[0]==-1/2

I strongly believe that its analytic solution can be written as an order-12 curve in $a$ and $u^2$ variables (implicit solution) $$ \sum_{i=0}^{6}\sum_{j=0}^{12}C_{i,j}u^{2i}a^j(u)=0, \text{where }C_{0,0}\neq0. $$ I know that Maple can solve such equations in principle (see comments below this post). But it is also known that even to verify such an implicit solution is difficult with MA (Maple seems to do the verification without problems).

I would very much love to obtain the solution with MA. I am expecting that most of the $C_{i,j}$ coefficients should be zero, and the remaining ones should be small positive or negative integer numbers. For this reason I used FindIntegerNullVector, however, didn't find anything meaningful. I suppose that numerical solution was not accurate enough.

Background

The differential equation above is the result of many transformations. The initial system of equations is presented in MO post is given here for convenience in Mathematica form

λ[a_]:=a^3+1/4*a
Λ[a_]:=1/4 a^4+1/8*a^2
sys3odes={2a'[u]==u*x'[u],
4u*a[u]*x'[u]-4y'[u]==-2x[u]*D[u*a[u],u],
4u*a[u]*y'[u]-4x'[u]==D[u*λ[a[u]],u]-2y[u]*D[u*a[u],u]};

By combining Eq.(I) and (II) and integrating with given initial conditions we obtain

eqy=2y[u]==u*a[u]*x[u]+a[u]^2-1/4;

By adding y * Eq.(II) and x * Eq.(III) and integrating with given initial conditions we obtain

eqxy=4 u*a[u]*x[u]*y[u]-2(x[u]^2+y[u]^2)-u*x[u]*λ[a[u]]+2Λ[a[u]]==2Λ[-1/2];

We can verify these algebraic relations numerically:

{sa, sx, sy} = NDSolveValue[
  Join[sys3odes, {a[0] == -1/2, x[0] == 0, y[0] == 0}], {a, x, y},
  {u, 0, 1000}, WorkingPrecision -> 100, MaxSteps -> Infinity]
Plot[-2 sy[u] + u*sa[u]*sx[u] + sa[u]^2 - 1/4, {u, 0, 10}]
Plot[4 u*sa[u]*sx[u]*sy[u] - 2*(sx[u]^2 + sy[u]^2) - u*sx[u]*λ[sa[u]] + 2 Λ[sa[u]] - 2 Λ[-1/2], {u, 0, 10}]

I do not present plots, but they are indeed very close to zero. Next you eliminate $x(u)$ and $y(u)$.

Eliminate[{sys3odes[[1]],eqy,D[eqy,u],eqxy,D[eqxy,u]},
{x[u],x'[u],y[u],y'[u]}]//SubtractSides//Factor
Collect[Expand[%[[1,3]],a'[u]],a'[u]]
CoefficientList[%,a'[u]]//Simplify

In this way one gets the $r,q,p$ coefficients presented at the top of my post.

Answer 1

Your equation does not have a unique solution since it is quadratic in $a'(u)$: it has two solutions.

One solution is $$ a(u)=-\frac{\sqrt{u^2-\sqrt{u^4-4 u^2+16}+4}}{\sqrt{6} u} $$ which satisfies the equation $$ -27 u^{10} a(u)^{12}+9 u^{10} a(u)^{10}+108 u^8 a(u)^{10}-33 u^8 a(u)^8-144 u^6 a(u)^8+40 u^6 a(u)^6+64 u^4 a(u)^6-16 u^4 a(u)^4=0 $$

Check:

eq /. a -> (-(Sqrt[4 + #^2 - Sqrt[16 - 4 #^2 + #^4]]/(Sqrt[6] #)) &) // Simplify
Limit[a[u] /. a -> (-(Sqrt[4 + #^2 - Sqrt[16 - 4 #^2 + #^4]]/(Sqrt[6] #)) &), u -> 0, Direction -> "FromAbove"]

(* {True, True} *)

I am not sure how to find the other solution, but I assume/hope this is enough for your purposes.

More generally, given arbitrary $p,q,r$, the curve you are after is $$ q(a,u)^2-4 p(a,u) r(a,u)=0 $$

Answer 2

Not a solution, but a start:

p[a_, u_] = -(16 + u^2)^2 + 3 u^2 a^2 (36 u^4*a^4 + 3*u^2*(u^2 - 48)*a^2 - u^4 + 8*u^2 + 192);
q[a_, u_] = 2*u*a*(9*u^4*a^4 - 3*u^2*(u^2 - 4)*a^2 - (u^2 + 16));
r[a_, u_] = u^2*a^2*(3*u^2*a^2*(3*a^2 - 1) - 1);

We can define the solution's series-expansion around 0, for example to order 200 in $u$:

Q = DifferenceRoot[Function[{y, n}, {2 n (1 + 2 n) (7 + 4 n) y[n] + (-135 - 308 n - 240 n^2 - 64 n^3) y[n+1] + (960 + 2144 n + 1344 n^2 + 256 n^3) y[n+2] == 0, 
                                     y[1] == 1/64, y[2] == 9/4096}]];
a[u_] = -1/2 + Sum[Q[n] u^(2n), {n, 200}];

This function satisfies the boundary condition,

a[0]
(*    -1/2    *)

and the differential equation to order 405 in $u$:

Series[r[a[u], u] + q[a[u], u]*a'[u] + p[a[u], u]*a'[u]^2, {u, 0, 405}]
(*    O[u]^406    *)

By studying the coefficients of the series expansion, you may learn more.

Array[Q, 5]
(*    {1/64, 9/4096, 17/131072, -491/16777216, -5049/536870912}    *)

For $n\gg1$ these coefficients are about $-0.115 \times 4^{-n} n^{-3/2} \cos[(n + 1) \pi/3]$:

DiscretePlot[Q[n]/(-4^-n n^(-3/2) Cos[(n + 1) π/3]), {n, 500}]

which means that the radius of convergence is only 2:

Plot[a[u], {u, 0, 2}]

How did I find the coefficients $Q$? I set a[u_] = -1/2 + Sum[q[i] u^(2 i), {i, 200}] (where the q[i] are unknown) and inserted this form into the differential equation, then series-expanded the differential equation around $u=0$ to high order, set all expansion orders to zero, and solved for the coefficients q[i]. Once I had enough of these, I used FindSequenceFunction to find a closed form.

Answer 3

Here an alternative solution approach

Clear[p, q, r]
eq = r[a[u], u] + q[a[u], u]*a'[u] + p[a[u], u]*a'[u]^2 == 0;
sol=Solve[eq, a'[u]]     
(*{{Derivative[1][a][u] -> (-q[a[u], u] - Sqrt[q[a[u], u]^2 - 4 p[a[u], u]r[a[u], u]])/(2 p[a[u], u])}, 
{Derivative[1][a][u] -> (-q[a[u], u] + Sqrt[q[a[u], u]^2 - 4 p[a[u], u] r[a[u], u]])/(2 p[a[u], u])}}*)

We only get a unique solution a[u] if the Sqrt-part q[a[u], u]^2 - 4 p[a[u], u] r[a[u], u] in sol vanishs!

null=q[a[u], u]^2 - 4 p[a[u], u] r[a[u], u]

With the definitions

p[a_, u_] = -(16 + u^2)^2 +3 u^2 a^2 (36 u^4*a^4 + 3*u^2*(u^2 - 48)*a^2 -u^4 + 8*u^2 + 192);
q[a_, u_] = 2*u*a*(9*u^4*a^4 - 3*u^2*(u^2 - 4)*a^2 - (u^2 + 16));
r[a_, u_] = u^2*a^2*(3*u^2*a^2*(3*a^2 - 1) - 1);

we get a list of possible solutions

sol = Solve[q[a[u], u]^2 - r[a[u], u] 4 p[a[u], u] == 0, a[u]]
(*{{a[u] -> 0}, {a[u] -> 0}, {a[u] -> 0}, {a[u] ->0}, 
{a[u] -> -(2/(Sqrt[3] u))}, {a[u] -> -(2/(Sqrt[3] u))}, 
{a[u] -> 2/(Sqrt[3] u)}, {a[u] -> 2/(Sqrt[3] u)}, 
{a[u] -> -(Sqrt[1 + 4/u^2 - Sqrt[16 - 4 u^2 + u^4]/u^2]/Sqrt[6])}, {a[u] -> Sqrt[1 + 4/u^2 - Sqrt[16 - 4 u^2 + u^4]/u^2]/Sqrt[6]}, 
{a[u] -> -(Sqrt[1 + 4/u^2 + Sqrt[16 - 4 u^2 + u^4]/u^2]/Sqrt[6])}, 
{a[u] -> Sqrt[1 + 4/u^2 + Sqrt[16 - 4 u^2 + u^4]/u^2]/Sqrt[6]}}*)

Only one of these fullfill the required condition ic

Limit[a[u] /. sol, u -> 0]
(*{0, 0, 0, 0, Indeterminate, Indeterminate, Indeterminate, \Indeterminate, -(1/2), 1/2, -\[Infinity], \[Infinity]}*)

pos=Position[Limit[a[u] /. sol, u -> 0], -1/2][[1, 1]]
sol[[pos]]

(*{a[u] -> -(Sqrt[1 + 4/u^2 - Sqrt[16 - 4 u^2 + u^4]/u^2]/Sqrt[6])}*)

$\left\{a(u)\to -\frac{\sqrt{-\frac{\sqrt{u^4-4u^2+16}}{u^2}+\frac{4}{u^2}+1}}{\sqrt{6}}\right\}$

That's the solution given in @AccidentalFourierTransform's answer unfortunately without derivation.