Series solution does not match numerical solution

I have a system of algebraic-differential equations (for background see this question): $$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.$$

with initial conditions $$a(0)=-1/2, x(0)=0$$. It can be nicely solved numerically

{ka, kx} =
NDSolveValue[{2 Derivative[1][a][u] == u Derivative[1][x][u],
4 a[u]^2 (1 + 3 u^2 x[u]^2) == 1 + 4 x[u] (u a[u] + 4 x[u]),
a[0] == -1/2, x[0] == 0}, {a, x}, {u, 0, 1000}];


Since I am interested in the behaviour of $$a(u)$$ at $$u\rightarrow\infty$$, and I suspect that the leading term of $$a(u)$$ must be $$-1/u$$, I seek a solution in the forms of series:

m = 7;
a[s_] := -1/s + Sum[A[i]/s^i, {i, 3, m, 2}]
x[s_] := Sum[B[i]/s^i, {i, 0, m, 2}]
eqs = Series[{2 Derivative[1][a][u] == u Derivative[1][x][u],
4 a[u]^2 (1 + 3 u^2 x[u]^2) == 1 + 4 x[u] (u a[u] + 4 x[u])}, {u, Infinity, m - 1}] // Normal;
{sol} = SolveAlways[eqs, u];


Now, the solutions can be compared:

Plot[{a[s] /. sol, ka[s]}, {s, 0, 20}]


and

Plot[{x[s] /. sol, kx[s]}, {s, 0, 20}]


. . . .

As you can see, the series solution for $$a$$ nicely match its numerical result. However, the series solution for $$x$$--does not match the respective numerical result. I would like to know if this is a problem with the default solver of NDSolveValue, or am I doing something wrong mathematically?

• Can you please re-check your code with NDSolveValue? I get a $x\equiv0, a\equiv-1/2$ solution ... Jan 11 at 16:17
• I assume you might have not cleared the definitions for x and a before running NDSolveValue. Jan 11 at 16:32
• @Domen I have checked. I really get what is in OP. My version is 13.1.0 for Mac OS X ARM (64-bit) (June 16, 2022). Jan 11 at 16:44
• I get an ndsz error after the first few steps. If I differentiate the algebraic equation, I get results that look like above. Changing precision and integration method does not alter the NDSolve results significantly. However, FWIW, changing the initial value of x[0] changes the limit that the numerical x[s] approaches: Plot[Evaluate@Table[NDSolveValue[{2 Derivative[1][a][u] == u Derivative[1][x][u], D[4 a[u]^2 (1 + 3 u^2 x[u]^2) == 1 + 4 x[u] (u a[u] + 4 x[u]), u], a[0] == -1/2, x[0] == x0}, x[s], {u, 0, 100}, Method -> Automatic], {x0, Range[0, 10]/10}], {s, 0, 20}] Jan 11 at 17:01
• I am on version 13.3.0 for Win64. I get the same result as Domen for the NDSolveValue (ie constants), and I can only reproduce the blue curve in the series result. Jan 11 at 17:04

[Update: I needed to recalculate the initial $$a(0)$$ for each different initial $$x(0)$$, so that the solution starts on the algebraic equation.]

If I rework the OP's approach with a differentiated algebraic equation, then the asymptotic solution is different and matches the numerical solution better. I don't understand the math. This sort of analysis is missing in my background. The following shows how well the asymptotics (color) match the numerics (dashed black), and (update) how well the algebraic equation is solved out to "infinity" (10^6):

I hope this is not cheating: The match shows that there is a numerical solution matching each asymptotic solution and vice versa. What happens is this. The asymptotics of differentiated system sol2 has a free parameter B[0], the limit at infinity. I solved for that parameter 11 times using the numerical solution for $$x(0) = 0.0, 0.1,\dots, 0.9, 1.0$$, and (update) solved for $$a(0)$$ using the algebraic equation. If the asymptotics are a good model, they should match, which they seem to. I don't know how to solve for the limit at infinity in terms of $$x(0)$$, except by numerically integrating a long way. (Update: Use of DSolveChangeVariables to switch $$u=\infty$$ to $$0$$ omitted, because of some silly activation issue means at the moment no access to V14/13.3. :/)

Sorry I can't explain the math better, but I hope this helps in some way. The comparison made it seem I was probably on the right track, and I wanted to share it in case it's helpful. Method -> {"Projection", "Invariants" -> {aeq}} in NDSolve does better at satisfying the algebraic equation, but there's no significant qualitative difference.

Code:

m = 7;
aa[s_] := -1/s + Sum[A[i]/s^i, {i, 3, m, 2}]
xx[s_] := Sum[B[i]/s^i, {i, 0, m, 2}]
eqs = Series[{
2 Derivative[1][a][u] == u Derivative[1][x][u],
D[4 a[u]^2 (1 + 3 u^2 x[u]^2) == 1 + 4 x[u] (u a[u] + 4 x[u]),
u]} /.
{a -> aa, x -> xx},
{u, Infinity, m - 1}] // Normal;
sol2 = SolveAlways[eqs, u];

xx[1000.] /. sol2 (* inspect solution for large u *)
(*
{-9.99997*10^-7 + 1.66667*10^-18 A[5],
-1.*10^-6 + 1. B[0] + 3.*10^-12/(1 + 6 B[0]) + (2.*10^-17 (1 + 7 B[0]))/(1 + 6 B[0])^3}
*)

ode = {2 Derivative[1][a][u] == u Derivative[1][x][u],
D[
aeq = 4 a[u]^2 (1 + 3 u^2 x[u]^2) == 1 + 4 x[u] (u a[u] + 4 x[u]),
u]};
inf = 10^6;
sols = Join @@ Table[NDSolve[{
ode,
{a[0] == (a[0] /.
FindRoot[aeq /. u -> 0 /. x[0] -> x0, {a[0], -1/2},
WorkingPrecision -> 40]), x[0] == x0}
}, {a, x}, {u, 0, inf}, WorkingPrecision -> 40,
InterpolationOrder -> All], {x0, 0, 1, 1/10}]; // AbsoluteTiming
(*{1.33258, Null}*)

Plot[{xx[s] /. sol2[[2]] /. B[0] -> (x[inf] /. sols), x[s] /. sols} //
Flatten // Evaluate
, {s, 0, 20}
, PlotStyle -> Join[
Table[Hue[3 i/(4 Length[sols])], {i, 0, Length[sols] - 1}],
Table[
Directive[AbsoluteThickness[1]; Dashed, Black], {Length@sols}]
]
, PlotLegends -> Append[N[x[0] /. sols], HoldForm[x]],
ImageSize -> 360, PlotLabel -> "Asymptotic vs. Numerical"]

LogLogPlot[
((xx[s] /. sol2[[2]] /. B[0] -> (x[inf] /. sols)) - (x[s] /.
sols))/(x[s] /. sols) // Abs // Evaluate
, {s, 0, inf/10}
, PlotStyle -> Join[
Table[Hue[3 i/(4 Length[sols])], {i, 0, Length[sols] - 1}],
Table[
Directive[AbsoluteThickness[1]; Dashed, Black], {Length@sols}]]
, PlotLegends -> N[x[0] /. sols], ImageSize -> 360,
PlotLabel -> "Relative Error", WorkingPrecision -> 40]

LogLogPlot[
(Subtract @@ aeq /. sols) // Abs // Evaluate
, {u, 0, inf}
, PlotStyle -> Join[
Table[Hue[3 i/(4 Length[sols])], {i, 0, Length[sols] - 1}],
Table[
Directive[AbsoluteThickness[1]; Dashed, Black], {Length@sols}]]
, PlotLegends -> N[x[0] /. sols], ImageSize -> 360,
PlotLabel -> "Residual Error in Algebraic Equation",
WorkingPrecision -> 40]

• I seems to understand your approach better. By solving two ODE you are not able to determine B[0] and therefore resort to numerics. My question would be, at $u=\infty$, does it fulfil the algebraic equation? Jan 12 at 12:14
• @yarchik A good point. You're probably right to doubt that. Jan 12 at 16:50
• @yarchik See if the update satisfactorily addresses your point. Jan 12 at 18:39
• Thanks for your update. There is a lot of information, but I think the conclusion should be that my series solution uses unjustified assumption that $a\sim-1/u$. Jan 13 at 7:38