Fourth-order BVP problem with boundary conditions at infinity

Question

I am trying to compute the solution of the fourth-order ODE

$$ y^{\prime \prime} + V(y) - \beta y^{\text{ (IV)} } = 0$$ with $V(x) = -2x + 4 x^3$, on the real line, with boundary conditions $$ \begin{cases} \lim_{x \pm \infty} y(x) = 0 \\ \lim_{x \pm \infty} y^ \prime (x) = 0 \end{cases}$$.

Clearly the trivial solution $ y = 0$ is there, but another solution should exist, with $y(0) \neq 0$. Exponential decay at infinity is expected, so the conditions $$ \begin{cases} \lim_{x \pm \infty} y ^ {\prime \prime} (x) = 0 \\ \lim_{x \pm \infty} y^ {\prime \prime \prime} (x) = 0 \end{cases}$$ should apply. On symmetry grounds, $y^ \prime (0) = 0$ is also expected.

I tried to solve it on an interval $(0,L)$ with NDSolve such as in

  L = 10; s = NDSolve[{y''[x] -2*y[x] + 4*y[x]^3 - y''''[x] == 0, y[L] == 0, 
  y'[L] == 0, y'[0] == 0, y''[L] == 0}, y, {x, 0, 10}, WorkingPrecision -> 60]

also trying other b.c., such as the third derivative, but no success

      FindRoot::cvmit: Failed to converge to the requested accuracy or precision within 100 
      iterations.
      NDSolve::berr: The scaled boundary value residual error of 2.9040010780303314`*^7 indicates 
      that the boundary values are not satisfied to specified tolerances. Returning the best solution 

I tried the Method -> {"Shooting", "StartingInitialConditions" -> {y`[0] == 0}}, using the boundary conditions at the extremity $L$ for NDSolve, to no avail

  NDSolve::mxst: Maximum number of 10000 steps reached at the point x == 
  6.2422167664683852553263167777483585328023492585909723902636`30..

I have searched the site and found similar "separatrix" problems addressed in the answers Computing separatrix of second order nonlinear autonomous ode, as well as Numerical solution of nonlinear boundary value problem, but I am unsure they are applicable. The former seems to me to rely on a particular structure for a system of second order ODEs, an the latter requires knowledge of the $y(0)$ value, if I am not missing the point.

Any hint would be most appreciated, thanks

EDIT

The differential equation has a physical interpretation as 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 on an elastic foundation $U(x) = x^2 - x^4$. The first term represents the stretching energy, the second the potential energy, the third the bending energy. I thought about using this fact, by adding a "guess" boundary condition $ y(0) = y_0$ and then optimise $y_0$ by making the energy, as defined by the functional above, stationary. But even choosing an arbitrary $y_0$ I struggle to solve the ODE, for example choosing two b.c. per end, function value and first derivative,

   s = NDSolve[{y''[x] - 2*y[x] + 4*y[x]^3 - y''''[x] == 0, y[L] == 0, 
   y'[L] == 0, y'[0] == 0, y[0] == 1}, y, {x, 0, 10}, 
   WorkingPrecision -> 30, MaxSteps -> 10000]

returns

  FindRoot::cvmit: Failed to converge to the requested accuracy or precision 
  within 100 iterations.
  NDSolve::berr: The scaled boundary value residual error of 
  3.0598614681511925`*^14 indicates that the boundary values are not 
   satisfied to specified tolerances. Returning the best solution found.

which I am puzzled about, I am not so sure what makes this BVP so challenging.

EDIT 2

I am trying to use the approach of heat convection differential equations from 1952 - Mathematica fails to converge. If it could be of any utility, the (interesting) solution of the above BVP for $\beta = 0$ is $$ y(x) = \operatorname{sech} (\sqrt{2}x)$$ as Mathematica confirms

           eq = u''[x] - 2 u[x] + 4 u[x]^3 == 0
           FullSimplify[eq /. u -> Function[{x}, Sech[ Sqrt[2] x]]]
           (*true*)

Answer 1

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

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

as follows

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

y0 = Exp[-Abs[x]];
beta = -0.00001;
xmax = 5;
nmax = 20;
error = 0.00001;
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]];
   y2 = y1[x] /. sol;
   solx = NMaximize[{Abs[y0 - y2], x >= -xmax, x <= xmax}, x][[1]];
   Print[{k, solx}];
   If[k == nmax || solx < error, thickness = Thick; color = Black; y0 = y2];
   AppendTo[SOLS, Plot[y0, {x, -xmax, xmax}, PlotStyle -> {thickness, color}]];
   If[solx < error, Break[], y0 = y2];
]

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]

Approaching $\beta$ by $-0$ with $\beta = -0.00001$ we obtain

Answer 2

If we look at the eigenvalues of the Jacobian at each value of y and beta, we see there is at least one eigenvalue that has a positive part and at least one that has a negative real part. So no matter which way you integrate, the numerical integration will be unstable. In fact most, if not all, solutions seem to develop singularities (poles). The tendency to blow up when integrating over long intervals is what makes solving the BVP so difficult. I cannot prove that a solution to the BVP does not exist, but it probably cannot be found numerically. I wonder if there is something missing or wrong in the ODE.

fop = Internal`ProcessEquations`FirstOrderize[
   {-2 y[x] + 4 y[x]^3 + y''[x] - \[Beta] D[y[x], {x, 4}] == 0}, {x}, 
   1, {y}];
Column[fop, Dividers -> All]
sys1o = Take[fop, 2] // Flatten; (* first order system *)
dvars = 
 Through[Flatten[fop[[3]]][x]];(* dependent variables *)
jac = D[ (* Jacobian *)
  First@Values@Solve[sys1o, D[dvars, x]],
  {dvars}
  ]
evals = Eigenvalues[jac] /. y[x] -> y (* local EVs as func. of y and beta *)

Manipulate[
 Block[{\[Beta] = b},
  ReImPlot[evals, {y, -0.8, 0.8}, MaxRecursion -> 10, 
   PlotStyle -> Table[AbsoluteThickness[4 - 2 k/3], {k, 4}], 
   PlotRange -> 100000, PlotRangePadding -> Scaled[.05], 
   ScalingFunctions -> {ArcTan, Tan}, Frame -> True,
   FrameTicks -> {{{-1000, -5, -1, 0, 1, 5, 1000}, Automatic}, 
     Automatic},
   FrameLabel -> {y, HoldForm@ReIm@\[Lambda]}]
  ],
 Row[{
   Control@{{b, 1/8, \[Beta]}, -17, 4, 
     Manipulator[Dynamic[Log2@b, (b = 2^#) &], #2] &},
   " ", Dynamic@Style[b, "Label"]}]
 ]

The equivalent first-order system is $$ \eqalign{ y_0' &= y_1 \cr y_1' &= y_2 \cr y_2' &= y_3 \cr y_3' &= L(y_0,y_1,y_2,y_3) + 4y_0^3/\beta \cr } $$ where $L(y_0,y_1,y_2,y_3)=-2y_0/\beta+y_2/\beta$ is the linear part of the ODE. When $y_0\approx0$, the system is approximately linear with a small perturbing term $4y_0^3/\beta$. At a given value of $y_0$, let $v_j(x)$ be the eigenfunctions of the linear part of the system with corresponding eigenvalues $\lambda_j$. If we write ${\bf y}=(y_0,y_1,y_2,y_3)^T$ and express the solution at a given value of $y_0$ by $${\bf y}(x) = \sum_j \alpha_j v_j(x) \,,$$ then we have $$d{\bf y} = \sum_j (\lambda_j\alpha_j+\delta_j) v_j(x) \,,$$ where the $\delta_j$ come from the perturbing term. It turns out that $\delta_j \ne 0$ (for $y\ne0$) although they are small. The terms associated with eigenvalues with $\mathop{\text{Re}}(\lambda_j)>0$ grow as the integration progresses. Now, it is possible that, since the value of $y_0$ changes and therefore the linear system and its eigenfunctions change, that the new $\delta_j$ might cancel the old one. I don't think so, not if $y_0 \rightarrow 0$, but I haven't proven it.

Answer 3

One possible way to solve this problem. Let put y[0]=y0 where $y0\ne 0$ is some parameter. Now we can normalize y->y0 yn so that yn[0]=1 and equation turns to this one $$yn''[x] -2*yn[x] + 4*y0^2 yn[x]^3 - yn''''[x] == 0$$.
This equation with parameter can be solve with ParametricNDSolveValue[] as follows (we omit n}

L = 10; s = 
 ParametricNDSolveValue[{y''[x] - 2*y[x] + 4*y0^2  y[x]^3 - 
     y''''[x] == 0, y[0] == 1, y'[0] == 0, y''[L] == 0, y'[L] == 0}, 
  y, {x, 0, L}, {y0}]

Finally we plot solution for a range of y0

Plot[Evaluate[Table[y0 s[y0][x], {y0, .05, .15, .02}]], {x, 0, L}, 
 PlotRange -> All, PlotLegends -> Table[y0, {y0, .05, .15, .02}]]