# Asymptotic solution of a system of ODEs

I have the following system of Ordinary Differential Equations (ODEs) together with initial values

ode = {-4 * (1 + u * a[u]) * f'[u] == (a[u]^3 + 1/4 * a[u] + (3 a[u]^2 + 1/4) * u * a'[u]) + 2 (a[u] + u * a'[u]) * f[u],
a'[u] == 1/2 * u * f'[u],
a[0] == -1/2, f[0] == 0}


which can be easily solved numerically:

{ra, rf} = NDSolveValue[ode, {a, f}, {u, 0, 20}]


and asymptotically for $$u\rightarrow0$$.

AsymptoticDSolveValue[odetph, {a[u], f[u]}, {u, 0, 3}]

(* {-(1/2) + u^2/64 + u^3/128, u/16 + (3 u^2)/128 + (3 u^3)/512} *)


But I need to know an asymptotic expansion of a[u] as $$u\rightarrow\infty$$. Even the coefficient of the lowest order term $$\sim 1/u$$ would be very nice to know. Is there a trick to obtain such a solution with MA?

A simple approach is to introduce the ansatz $$f(u)=\sum_{k\ge0}\frac{f_k}{u^k},\qquad a(u)=\sum_{k\ge0}\frac{a_k}{u^k}$$ which you plug into the ODE and solve for the coefficients. I find \begin{aligned} f(u)&=\frac{3}{8}-\frac{1}{u^2}+\frac{3}{2 u^4}-\frac{7}{2 u^8}+\frac{33}{2 u^{12}}+\cdots\\ a(u)&=-\frac{1}{u}+\frac{1}{u^3}-\frac{2}{u^7}+\frac{9}{u^{11}}+\cdots \end{aligned}

Check:

ode[[1 ;; 2]] /. {f -> (3/8 - 1/#^2 + 3/(2 #^4) - 7/(2 #^8) + 33/(2 #^12) &), a -> (-(1/#) + 1/#^3 - 2/#^7 + 9/#^11 &)}
Series[%, {u, ∞, 12}] // Simplify

(* {True, True} *)


You can of course add more terms, with symbolic coefficients, and find their correct value by solving the equation order by order. I'm sure you know what I mean but feel free to ask, of course.

The numeric minus the asymptotic solutions differ by

which is so close to 0 that the numerics overwhelm the error. Might need to increase accuracy to get a meaningful plot...

• Amazing! I was so focused on AsymptoticDSolve and why it does not work that I started thinking there is some fundamental obstacle and even did not try to use such an ansatz. Do you know why AsymptoticDSolve does not work here? Also, I understood you obtain the solution "by hands". Is there a way to automatise it? Jan 2 at 17:13
• @yarchik: You asked " Is there a way to automatise it?". Maybe, SolveAlways helps with it. Jan 2 at 17:26
• @user64494 I tried SolveAlways, but not all coefficients can be obtained in this way. I am not getting $f_0=3/8$ and $a_1=-1$ in this way. Jan 3 at 8:11
• @yarchik: ode = {-4*(1 + u*a[u])* f'[u] == (a[u]^3 + 1/4*a[u] + (3 a[u]^2 + 1/4)*u*a'[u]) + 2 (a[u] + u*a'[u])*f[u], a'[u] == 1/2*u*f'[u]} /. {a[u] -> Sum[a[k]/u^k, {k, 0, 3}], f[u] -> Sum[f[k]/u^k, {k, 0, 3}], a'[u] -> D[Sum[a[k]/u^k, {k, 0, 3}], u], f'[u] -> D[Sum[f[k]/u^k, {k, 0, 3}], u]};SolveAlways[ode, u] works for me. Jan 3 at 10:42
• @AccidentalFourierTransform Thank you. In the meantime I was able to solve the equations analytically. The search series can then be obtained as AsymptoticSolve[a^2 (2+a u)^2==1+a u,a,{u,Infinity,12}] Jan 3 at 21:11