# Analytic solution expected

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.

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$$

• I have to say, this is the most amazing answer on StackExchange that I have ever seen! I encourage you to provide a bit more details about your last equation and to answer and collect bonus at Mathoverflow. Jan 17 at 15:16
• @yarchik Thank you for the nice words! I'm not sure I understand the relation between this post and the one on overflow. If you take the solution here and expand around $u=\infty$ you find $a\sim -1/u+\cdots$ which is different from the asymptotics $a\sim-1.024/u+\cdots$. The two systems do not seem to be equivalent, no? Jan 17 at 15:49
• See my derivations in OP. Jan 17 at 16:28
• @yarchik Your initial ODE has a unique solution but your manipulations introduced a second, fictitious solution (as I mentioned above, your reduced system has two solutions). I only wrote one of the two solutions and I suspect that I chose the wrong one. What I wrote above is most likely the fictitious solution, not the one you want. I need to think whether it is possible to find the other solution as well. Jan 17 at 16:42
• @AccidentalFourierTransform : I don't understand your last sentence: "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$." If e.g. $p(a,u)=1$, $q(a,u)=-2$, and $r(a,u)=1-u^2$ (for $u\ge0$), then $a(u)=c+u\pm u^2/2$ for some constant $c$, whereas "the curve" $q(a,u)^2-4 p(a,u) r(a,u)=0$ is $u=0$. Jan 19 at 17:01

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.

• You said this can be a start. What would be the next step, what is the method to get an analytic solution (most probably implicit) from the series expansion? Maybe you can illustrate your idea with a simpler equation, such as in the linked post? Jan 16 at 23:07
• Next step would be to take my solution $a(u)$ and insert it into your proposed implicit form, then solve for the coefficients $C_{i,j}$. I tried it and it does not work: I don't think your proposed implicit solution can be fulfilled. Jan 17 at 9:23

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.

• Thank you very much. Actually this is a question to both of you (@AccidentalFourierTransform), I have never seen anywhere that a solution of an ODE can be obtained by demanding that derivative is unique. Is it well known? Does it follows from the uniqueness of the solution of initial value problem (Picard–Lindelöf theorem)? Jan 17 at 16:33
• @yarchik I (==engineer) thought much more simply and used a unique differential equation as a basis to get a unique solution. Jan 17 at 18:46