Asymptotic solution to the differential equation $y”+(\pi^2+\varepsilon)y+y^2=0$ with $y'(0)=y'(1)=0$

Question

In studying perturbation methods, I am stuck with the following problem.

Let $\varepsilon>0$ be a small parameter. Find the asymptotic approximation to the non-constant solution of the differential equation
$$
y''+(\pi^2+\varepsilon)y+y^2=0,\\
y'(0)=y'(1)=0.
$$

I tried to solving it by perturbation method. Assume that we have the expansion
$$
y(x)=\varepsilon^\alpha y_1(x)+\varepsilon^\beta y_2(x)+\cdots, \quad (*)$$
where $0<\alpha<\beta$. Substituting it into the original differential equation, at leading order $O(\varepsilon^\alpha)$, we obtain
$$
y''_1+\pi^2 y_1=0,\\
y'_1(0)=y'_1(1)=0,
$$
which yields $y_1(x)=A\cos(\pi x)$, where $A$ is a constant. To find the value of $A$, we need to look at the next order, which is
$$
\varepsilon^\beta y''_2+\varepsilon^\beta y_2+\varepsilon^{1+\alpha }y_1+\varepsilon^{2\alpha} y_1^2=0,\\
y'_2(0)=y'_2(0)=0$$
I tried to balance the terms by choosing $\beta=1+\alpha$ or $\beta=2\alpha$ or $1+\alpha=2\alpha$ or $\beta=1+\alpha=2\alpha$. However, in each of the cases, I ended up with $A=0$ when substituting the solution of $y_2$ into the boundary conditions.

On the other hand, integrating the original differential equation once, we have
$$
y'^2+(\pi^2+\varepsilon)y^2+\frac{2}{3}y^3=C,$$
where $C$ is a constant. Using the boundary conditions $y'(0)=y'(1)=0$, we can find the values of $y(0)$ and $y(1)$ once $C$ is specified. By separating variables, we see that
$$
\int_{y(0}^{y(1)}\frac{dy}{\sqrt{C-(\pi+\varepsilon^2)y^2-\frac{2}{3}y^3}}=\int_{0}^1 dx=1,$$
which can be used to find the constant $C$. Numerical calculations show that for small $\varepsilon>0$ , $C$ exists, so the differential equation has non-constant solutions for small $\varepsilon>0$. For instance, when $\varepsilon=0.01$, $C\approx 1.1699$ and the graph of a non-constant solution plotted by Mathematica is

I think that problem is that the asymptotic expansion of the non-constant solution may not be of form (*), but I don't know how what the correct form is.

Answer 1

Equivalent forms of the equation

With the substitution $y=Au(Bt)$ scaling both the domain and the range, the equation reduces to $$ u''+\frac{\pi^2+ε}{B^2}u+\frac{A}{B^2}u^2=0. $$ The scales can now be chosen in a way that produces the unperturbed equation with $ε=0$ using $A=B^2=1+\frac{ε}{\pi^2}$, or a form closer to the standard harmonic oscillator $$ x''+x+δx^2=0 $$ with $B^2=\pi^2+ε$, $A=δB^2$.

Citing a previously performed perturbation analysis

This last form, after switching the sign in $x$, was explored in a previous topic. Fixing the amplitude $R$, the solution close to the circle of radius $R$ is given by the parametrization $$ x\sqrt{1+\frac23δx}=-R\cos\phi, ~~ \dot x=R\sin\phi,\\ \implies x=-R\cos\phi-δ\frac{R^2}3\cos^2\phi+O(δ^2),\\ \left[\small\begin{aligned} \implies 1&=(1+\frac{2δR}3\cos\phi+...)\dot\phi,\\ t&=\phi+\frac{2δR}3\sin\phi+...\iff \phi=t-\frac{2δR}3\sin t+...\\ x(t)&=-R\left[\cos t+\sin t\left(\frac{2δR}3\sin t\right)+...\right]-δ\frac{R^2}3\cos^2t+O(δ^2)\\ &=-R\cos t-δ\frac{R^2}3[\cos^2t+2\sin^2 t]+O(δ^2) \end{aligned}\right] $$ for the half-period of interest $\phi\in[0,\pi]$. Taking the derivative of the first equation one can derive a first-order ODE for the angle. Doing some series expansion one gets an expression for the half-period $\phi(T)=\pi$ as $$ T=\pi\left(1+\frac{5}{12}(δR)^2+\frac{385}{576}(δR)^4+O(δ^6)\right). $$ Per the parametrization we want $B=T$, which in the leading order gives $$ (\piδR)^2=\frac65ε. $$ This suggests $δ=\sqrtε$, $R=\frac{\sqrt{6/5}}{\pi}$.

Fall back to perturbation series for numerics

Now set $x=-R\cos t+δx_1$ to result in the equation $$ x_1''+x_1+x^2=0, $$ which in first order has the initial condition $x_1(0)=-\frac{R^2}3$ and the solution $$ x_1=\frac{R^2}6(\cos(2t)-3). $$ This can be used in an approximation $x=x_0+δx_1+δ^2v$ or just as initial guess in the boundary value problem for $v$ in $x=x_0+δv$.

Numerical experiments

# global constants and invariant functions
R2 = 1.2/np.pi**2
R = R2**0.5

def x0(t): return -R*np.cos(t) 
def dx0(t): return R*np.sin(t)
def x1(t): return R2/6*(np.cos(2*t)-3)
def dx1(t): return -R2/3*np.sin(2*t)
    
# Variables depending on the perturbation parameter
eps = 1e-8
dta = eps**0.5
B2 = np.pi**2+eps
B, A = B2**0.5, dta*B2
dvb = -dx0(B)/dta
 
# Boundary Value Problem
def ode_v(t,v): xt= x0(t)+dta*v[0]; return [v[1], -v[0]-xt**2]
def bc_v(v0,vb): return [v0[1],vb[1]-dvb]
# initial guess
t = np.linspace(0,B,5)
v = [ x1(t), dx1(t) ]
# solve and print (error) message
res = solve_bvp(ode_v,bc_v,t,v, tol=1e-4)
print(res.message,res.y[0,0], res.y[0,-1],-R2/3)

This converges well and reliably. Plotting after scaling back for the original problem gives in the phase plane