# Fourth-order BVP problem with boundary conditions at infinity

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.242216766468385255326316777748358532802349258590972390263630..


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*)

• It looks like we need to put bc as y''[0]==1, and not as y''[L]==1. Is it correct? Jan 18, 2021 at 13:46
• @Alex Trounev, that is indeed a typo, I meant there y'' [L] = 0. I think nothing can be said for the second derivative at the origin, while the condition y' [0] = 0 is expected due to symmetry, on top of all the b.c. at infinity. Jan 18, 2021 at 13:52
• But we can't put only zeros on every border since $y=0$ is solution. Jan 18, 2021 at 13:57
• @Alex Trounev, that is exactly where the problem is! I do not know of any other, non-zero, boundary condition. I will edit the post to provide some physical background, maybe it could be useful. Jan 18, 2021 at 14:08
• If you are expecting an even solution ($y'(0)=0$) then applying $y'''(0)=0$ makes sense. Jan 18, 2021 at 14:13

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

• thanks for this, very instructive code. It solves the case $\beta = 1$ very quickly, while NDSolve seem to struggle, as I use it. The fact is, the value $y(0)$ is not known a priory, and cannot hence be used as a b.c. I am sorry if I failed to make this point clear in my OP. The b.c. I have lead inexorably to the trivial solution, but that is it, is there a way to find the interesting solution. Your code is very useful and could be used to optimise over $y(0)$, I am still wondering if there is anything better. I will start a bounty, failing which I will certainly accept our answer. Jan 20, 2021 at 20:38
• @Smerdjakov I modified the script to handle two sided boundary conditions and seems to work well. Jan 20, 2021 at 21:00
• that is actually great I will look into it in depth, gives plenty of insight, how the iteration proceeds to the non-trivial solution using your absolute value exponential, thanks Jan 20, 2021 at 21:03

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 = InternalProcessEquationsFirstOrderize[
{-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.

• many thanks for this, very interesting. Yes the solution is unstable and that was expected, was hoping that with some magic, good initial guess, working precision, clever algorithms, it could be done. I am just wonder what do you mean by "the ODE is wrong"? You mean its instability suggests it is likely not to represent a physical system? It comes as Euler-Lagrange eq. of a pretty standard energy functional. Jan 21, 2021 at 8:11
• @Smerdjakov The instability is not just a numerical artifact. It implies that for almost all (or for all) initial conditions, $y$ as well as its derivatives $y^{(k)}$ go to $\pm\infty$, which implies the energy goes to infinity. I would have posted this a few days ago, but I felt I might have made a mistake. All I can say is that it seems consistent with the numerical problems we have with the BVP and why it’s feasible to solve it for xmax = 5 but not for xmax = 50. It made me wonder about the ODE. Alternatively, if the ODE is good, I was hoping someone would point out where I went wrong. Jan 21, 2021 at 13:36
• thanks for your reply. Indeed the instability is physical, that was expected. There is one a particular solution that is an unstable fixed point, and which is physically very interesting. I simply did not understand how badly such physical instability would have shown in numerical attempts, on which I do not know that much. Thanks again very useful. Jan 21, 2021 at 13:55

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}]]


• many thanks for this, it is very interesting how the normalisation makes it work. Jan 21, 2021 at 15:07
• @Smerdjakov you are welcome! It is common regularization for nonlinear model with a small parameter. Jan 21, 2021 at 21:09