Non-unique solutions for ODE with boundary conditions at infinity

Question

Ello,

I am looking for solutions to equations such as

$$ U ^ \prime (y) – y^{\prime \prime} + \beta y^{\textit{IV}} = 0$$
where $\beta >0$ is a constant and $U(x)$ is a function whose Taylor expansion at the origin equals $ \alpha x^2 + \dots$ for $\alpha >0$ and such that $\lim_{x \to \infty} = -\infty$.
A prototypical example would be the quartic double-well, $$ U(x) = x^2 -x^4$$.

There is a trivial solution $y(x) = 0$. On intuitive grounds (please see physical background below), as well as noting that the origin is a hyperbolic point for a dynamical system governed by the related system of equations, another solution is expected as well, with boundary conditions $$ \begin{cases} \lim_{x \pm \infty} y(x) = 0 \\
\lim_{x \pm \infty} y^ \prime (x) = 0 \end{cases}$$

This fact does not contradict the unicity theorem, as far as I understand, as the latter deals with boundary conditions on the real line, I believe.

Now to my problem, what numerical scheme could be used to get the non-trivial solution? Anything I could think of will pick the trivial solution branch only.
Further, I would be grateful for theory references on the possibility of multiple solutions for boundary conditions at infinity.

Things would be different if $\beta$ were equal to $0$.
In that case it is easy to determine two boundary conditions (also clarified in the physical background section below) for $x = 0$, and any numerical method will handle it easily.

Physical Background

The ODE is the Euler-Lagrange equation for the functional

$$ \int _{- \infty} ^{\infty } \frac{1}{2} y^{\prime 2} + U(y) + \frac{\beta}{2} y^{\prime \prime 2} \mathrm{d}x $$
representing the energy of a stretchable beam in the potential $U$. The first term represents the stretching energy, the second the potential energy, the third the bending energy.
Given this physical picture, the trivial solution represents a beam lying in the potential energy "valley" at the origin, un-stretched and straight. The non-trivial solution represents a beam lying on the valley as before, then deforming over the potential energy "bump", kept in the equilibrium by the pull from the descending potential energy on the other side, and going back to the valley at the origin.

If $\beta = 0$ the non-trivial solution is amenable to a simple solution.
The energy functional reduces in such case to

$$ \int _{- \infty} ^{\infty } \frac{1}{2} y^{\prime 2} + U(y) \mathrm{d}x $$ which could be interpreted as the action of a particle rolling in a potential energy $-U(x)$ (with $x$ now representing time). The turning point is then easily found by solving the equation $ U(x_1) = 0$: that is, the particle starts from $x = 0$, rolls down and up until $x_1$, and rolls back in a homoclinic orbit.
This gives two boundary conditions

$$ \begin{cases} y(0) = x_1 \\
y^ \prime (0) = 0 \end{cases}$$

and these allow any numerical scheme to pick the non-trivial branch.

If $\beta \neq 0$, this simple reasoning fails for the full fourth-order ODE (I also thought about converting the fourth-order ODE to a system of equations, to attempt a similar interpretation for a system of four particles, but did not get anywhere yet, apart confirming that a second non-trivial solution exists). I am hence left with the boundary conditions a infinity, and cannot get any solution other than the trivial one.

To sum up, I would like to be pointed towards references for the non-unicity of solutions when boundary conditions at infinity are applied, and to understand how to numerically solve similar problems.

The ODE has a physical interpretation and I was unsure whether the question were more apt for Physics StackExchange, or even SciComp StackExchange. I would certainly re-route if advised, I thought to start here as ultimately the question is primarily on the maths, thanks a lot.

Answer 1

You can use a linearization procedure to solve this nonlinear boundary value. For instance using a linear approximation scheme consisting in iterations over the linear ODE

$$ y''_{k+1}-V(y_k)+V'(y_k)(y_{k+1}-y_k)-\beta y''''_{k+1}=0 $$

given $y_0$ and the boundary conditions

$$ \cases{ y_{k+1}(-x_{max})=0\\ y'_{k+1}(x_{max})=0\\ y''_{k+1}(-x_{max})=0\\ y''_{k+1}(x_{max})=0 } $$

Let us proceed with the calculations for $\beta = 0$ which is known to have as solution

$$ y=\sqrt{1-\tanh ^2\left(\sqrt{2} x\right)} $$

We will approach this solution using negative $\beta$'s as follows in the attached MATHEMATICA script.

Clear[V, dV, y1, y0]
V[x_] := -2 x + 4 x^3;
dV[x_] := -2 + 12 x^2

y0 = Exp[-x];
beta = -0.00001;
xmax = 5;
nmax = 20;
bcs = {y1[-xmax] == 0, y1[xmax] == 0, y1'[-xmax] == 0, y1'[xmax] == 0};
SOLS = {Plot[y0, {x, -xmax, xmax}]};
thickness = Thin;
color = Blue;

For[k = 1, k <= nmax, k++,
   ode = y1''[x] + V[y0] + dV[y0] (y1[x] - y0) - beta y1''''[x];
   sol = NDSolve[Join[{ode == 0}, bcs], y1, {x, -xmax, xmax}][[1]];
   y0 = y1[x] /. sol;
   If[k == nmax, thickness = Thick; color = Black];
   AppendTo[SOLS, Plot[y0, {x, -xmax, xmax}, PlotStyle -> {thickness, color}]]
]

grfin = SOLS[[Length[SOLS]]];
Show[Plot[Sqrt[1 - Tanh[Sqrt[2] x]^2], {x, -xmax, xmax}, PlotStyle -> {Thick, Red}], grfin]
Show[SOLS, grfin, PlotRange -> All]

so with $\beta = -0.00001$ we obtain