# Numerical bounce solution

I want to solve numerically the following differential equation

$$y''(x) + \frac{3}{x}y'(x) = \frac{d U(y)}{y}\,,\qquad U(y) = \frac{1}{4}y(x)^2(y(x)-1)^2 -\frac{1}{400}y(x)^3$$

in the region $$x\in[0,+\infty]$$ with the bounce boundary conditions

$$y'(0) = 0\,,\qquad y(+\infty) = 0.$$

This problem is taken from Coleman's paper (click here) where he claims that such solution always exists, see page 4, right-column independently of the expression of the potential $$U(y)$$.

I want to find such solution for the potential that I have chosen.

This is a boundary-values problem and it is better to turn it into a initial-values problem. Let us then estimate the size of $$y(0)$$. It will become clear that I expect a flat direction around $$x=0$$, so substituting $$y'(0)=0$$ and $$y''(0) = 0$$ in the differential equation, I get an estimate for the initial condition $$y(0) \sim 1$$. This is just an estimate. I don't want the exact solution to satisfy $$y''(0)=0$$ identically, but I expect the second derivative to be small.

We can try to guess the form of the solution by making some approximations. Indeed, I can solve analytically the differential equation for big values of $$x\geq R \gg 1$$, by neglecting the $$1/400y(x)^3$$ term in the expression of $$U(y)$$, as $$y(0) \sim 1$$ and $$y(+\infty) = 0$$; then, by defining $$U_0(y) = \frac{1}{4}y(x)^2(y(x)-1)^2$$ we can solve

$$y''(x)= \frac{d U_0(y)}{y}\,$$

from which we get

$$y_R(x) = \frac{1}{2}\left(1-\text{tanh}\left[\frac{x-R}{2\sqrt{2}}\right]\right).$$

So, roughly the approximated solution is

$$y(x) = \left\{\begin{array}{c c} 1 & x \ll R\\ y_R(x) & x \geq R \end{array}\right.$$ and it is shown in the following picture for $$R=400$$

So, I have tried to solve numerically this Cauchy problem without any approximation by setting a cut-off $$R$$ on the $$x$$ coordinates and using the following code

V[y_[x_]] = 1/4 (-1 + y[x])^2 y[x]^2 - 0.0025 y[x]^3

Psol = ParametricNDSolveValue[{D[y[x], {x, 2}] + 3/x D[y[x], x] == D[V[y[x]], y[x]], y'[10^-5] == 0, y[10^-5] == tp, WhenEvent[y[x] < 0 || y[x] > tp, "StopIntegration"]}, y, {x, 10^-5, R}, {tp, R}]

ListPlot[Table[{tp, Abs[Psol[tp, 50][50]]}, {tp, 0.8, 1.015, 0.001}]]


where the last line may be useful to determine the magnitude of tp such that I recover a vanishing asymptotic value. However, I get the following plot

So, I guess the method is not suitable for this problem. In particular, the solution Mathematica finds does not respect the asymptotic boundary condition $$y(+\infty) = 0$$ at $$x=50$$, and moreover it starts oscillating after some time like in the following picture and does not go to zero at infinity

How can I get the numerical solution?

• My recommendation for you is to solve for your bounce using the shooting method described in this answer. Jan 29, 2019 at 19:03
• @QuantumDot thanks. I have tried to use shooting methods but without any success. BTW I am not an expert of numerical methods. Jan 30, 2019 at 13:44
• The initial value problem solver of NDSolve is very robust, the result given by your code only suggests one thing in my view: there just exists no solution satisfying $y(0)∼1,y′(0)=0,y(+∞)=0$. Jan 30, 2019 at 16:28
• @xzczd thank you for the comment. I have also thought about the non-existence of the solution. However, the fact that the approximate solution does exist and the term $\sim 1/400 y(x)^3$ is very small, makes me hoping that the exact numerical solution does exist as well. Is my impression too naive? Jan 30, 2019 at 16:55
• @MichaelE2 This problem is taken from Coleman's paper journals.aps.org/prd/abstract/10.1103/PhysRevD.15.2929 where he claims that such solutions always exists, see page 4, right-column Jan 31, 2019 at 18:17

## 3 Answers

The argument in the paper linked in the comments suggests that if one starts close enough to the one equilbrium at $$y \approx c_1 = 1.04$$, then the asymptotically vanishing solution may be found. Below is a solution that is not quite there. It starts at $$x = 10^{-8}$$ at $$y = c_1 - 10^{-25}$$ and uses WorkingPrecision -> 50; it certainly needs more than machine precision to ensure that $$y$$ does not round to $$c_1$$. It's not close enough. Getting close enough seems to require a very small difference and consequently a very high working precisions. But NDSolve fails when I push it, indicating an error test failure (NDSolveValue::nderr). It may be that for such an example to be computationally feasible, you need a greater difference in the value of V[x] at the minima.

ClearAll[V];
V[y_] = 1/4 (-1 + y)^2 y^2 - 25/10000 y^3;
ode = y''[x] + (3 y'[x])/x == V'[y[x]];

wp = 50;
cps = NSolve[V'[y] == 0, WorkingPrecision -> wp];

ndsol[x0_?NumericQ, x1_?NumericQ, y0_?NumericQ, yp_?NumericQ] :=
NDSolveValue[{ode, y'[x1] == yp, y[x1] == y0,
WhenEvent[y[x] < 0, "StopIntegration"],
WhenEvent[y[x] > 49/100 && x > x1, "StopIntegration"]},
y, {x, x0, 10000},
WorkingPrecision -> Precision@{x0, x1, y0, yp}];
ListLinePlot[
ndsol[2^-50, 1/10^8, y - 10^-25 /. Last@cps, 0],
PlotRange -> {All, {-1, 2}}]


Update

Here is a better result, FWIW. The method "Extrapolation" is more robust on this problem.

wp = 150;  (* working precision *)
cps = NSolve[V'[y] == 0, WorkingPrecision -> wp];
iObj[x0_?NumericQ, x1_?NumericQ, y0_?NumericQ, yp_?NumericQ] :=
Catch@NDSolveValue[{ode, y'[x1] == yp, y[x1] == y0,
WhenEvent[y[x] < 0, Throw[x]],
WhenEvent[y[x] > 49/100 && x > x1, Throw[$$Failed]], WhenEvent[y[x] > 2 && x > x1, Throw[$$Failed]]}, y, {x, x0, 1000},
WorkingPrecision -> Precision@{x0, x1, y0, yp},
Method -> "Extrapolation"];
obj[dy_?NumericQ] := iObj[2^-50, 1/10^8, y - 10^dy /. Last@cps, 0];

foo = Nest[  (* bisection method for finding optimal IC *)
With[{mid = obj[Mean@#]},
Switch[mid,
$Failed, {First@#, Mean@#}, _?NumericQ, {Mean@#, Last@#}, _, {First@#, Mean@#} ] ] &, {-42 - 54/100, -42 - 53/100}, 60 ] (* {-(2452178951689226635963/57646075230342348800), -(196174316135138130877/4611686018427387904)} *) ndsol[x0_?NumericQ, x1_?NumericQ, y0_?NumericQ, yp_?NumericQ] := NDSolveValue[{ode, y'[x1] == yp, y[x1] == y0, WhenEvent[y[x] < 0, "StopIntegration"], WhenEvent[y[x] > 49/100 && x > x1, "StopIntegration"], WhenEvent[y[x] > 2 && x > x1, "StopIntegration"]}, y, {x, x0, 10000}, WorkingPrecision -> Precision@{x0, x1, y0, yp}, Method -> "Extrapolation"]; ListLinePlot[ ndsol[2^-50, 1/10^8, y - 10^dy /. Last@cps, 0] /. dy -> First@foo, PlotRange -> All]  • Ehi Michael E2, thank you! I am trying to reproduce your result. Which version of Mathematica are you using? I am not able to reproduce your plots, by just copy and past your code, see here the screenshot ibb.co/4jdM6Fv Feb 1, 2019 at 8:59 • @apt45 I think the problem is that a couple of days ago I changed the def. of V[] so that the derivative would work with prime. I forgot about that. I included the new def. of V[] just now. If that doesn't work, I'm using V11.3 on a Mac. I doubt that would be a problem. Feb 1, 2019 at 11:38 • Thanks! Now I am able to produce the first plot, whereas regarding the second I get Error test failure like in the screenshot ibb.co/9hzgnJR . It may be a problem with the Working precision? Feb 1, 2019 at 11:47 • @apt45 Sorry again. I posted just before going to bed. I changed ndsol[] to use "Extrapolation"; code included. The default method (LSODA) gives those conv. test errors. Feb 1, 2019 at 11:53 • @apt45 Two issues with the working precision are that (1) the initial condition$y_0=c_1-10^{-p}$should be a distinct number from the critical value$c_1$and (2) the difference$y_0-c_1$should have sufficient precision. For (1) the WP should be more than$p$and for (2) the rule of thumb is twice that or$2p$. Since I didn't know$p\$ at the beginning, I put a 1 in front of the 50 from the first answer to get wp = 150, which turned out to be good enough. Beyond that, because the asymptotic solution is unstable, an unknown amount extra precision is needed to keep the solution on track. Feb 2, 2019 at 16:52

(Not an answer, an extended comment.)

So, I guess the method is not suitable for this problem. How can I get the numerical solution??"

Those statement and question do not make much sense to me. You already have the numerical solution. Why do you think the method is not suitable? The message InterpolationFunction::dmval is avoided if you decrease the range of tp a little (with 0.001). Study this result plot:

ListPlot[Table[{tp, Abs[Psol[tp, 50][50]]}, {tp, 0.8, 1.014, 0.001}], PlotRange -> All]


Also, maybe you should study the solution behavior for different values of the independent variable x not just 50.

For example, the following 3D plot gives an idea of the found solution behavior with respect to both tp and x. (We can see the "bouncing" along the x-axis.)

points = Flatten[Table[{tp, x, Abs[Psol[tp, 50][x]]}, {tp, 0.8, 1.014, 0.001}, {x, 1, 50, 1}], 1];
ListPlot3D[points, AxesLabel -> {"tp", "x", "y"}, PlotRange -> All]


• I am requiring my function to decay at infinity (see the second boundary condition). The solution you claim to have found does not seem to decay to zero for large x, even if I choose R=600 Jan 30, 2019 at 14:56
• @apt45 "The solution you claim to have found" -- read carefully what I wrote. I am not claiming to have found any solution, I am using your solution found with the code you posted. Jan 30, 2019 at 14:58
• Indeed! It is not the solution with the boundary conditions I have imposed! This is why I said that my method was not suitable Jan 30, 2019 at 14:59
• BTW I know you have copied and past my code, and then the function is the same that I got. For this reason, the question is still open, as I am not able to find the solution to the Cauchy problem in the OP Jan 30, 2019 at 15:01
• Anton have you seen the updated question? Jan 31, 2019 at 21:35

OK, Michael E2's answer gives me the courage to post the following as an answer. It seems that finite difference method (FDM) can be used to find the solutions for this problem in a relatively easy way. I'll use pdetoae for the generation of finite difference equations:

Clear[V];
V[y_] = 1/4 (-1 + y)^2 y^2 - 25/10000 y^3;
ode = y''[x] + (3 y'[x])/x == V'[y[x]];
newode = # x & /@ ode // Simplify
inf = 50;
domain = {0, inf};
bc = {y'[0] == 0, y[inf] == 0};
points = 500;
difforder = 2;
grid = Array[# &, points, domain];
(* Definition of pdetoae isn't included in this post,
please find it in the link above. *)
ptoafunc = pdetoae[y[x], grid, difforder];
ae = ptoafunc[newode][[2 ;; -2]];
aebc = ptoafunc@bc;

initial[x_] := 2;
sollst = FindRoot[{ae, aebc}, Table[{y@x, initial@x}, {x, grid}]][[All, -1]];
ListLinePlot[sollst, DataRange -> domain]
sol = ListInterpolation[sollst, grid];
(* Error check: *)
Plot[Subtract @@ ode /. y -> sol // Evaluate, {x, 0, inf}, PlotRange -> All]


I hesitated to post this because:

1. Using different inf, I can find different solutions, e.g. inf = 100:

2. The error of the solution isn't that small, and even if I increase points, I cannot observe obvious error reduction:

But now, reading all the analysis above, I believe there just exist multiple solutions for the problem, and the code above is worth sharing.

• I think this is reasonable because theoretically we know there's a solution, so you know you're approximating something that really exists. OTOH, the sensitivity to the truncation error near x == inf (I assume) seems to affect the computed solution quite a lot. Feb 2, 2019 at 17:14