# An algorithm to obtain appropriate numerical solutions for the bounce

Posted 9 months ago
719 Views
|
3 Replies
|
1 Total Likes
|
 I am looking to replicate an algorithm explained in the paper "Impact of new physics on the EW vacuum stability in a curved spacetime background" by E. Bentivegna, V. Branchina, F. Continoa and D. Zappal to numerically solve the boundary-value problem encountered when trying to find the bounce solution in false vacuum decay. It can be found here: https://arxiv.org/abs/1708.01138In the Appendix, "A Numerical computation of the bounce solution", it outlines an algorithm involving a shooting method to solve the following boundary problem: $y^{\prime\prime}(x) + \frac{3}{x}y^{\prime}(x) = \frac{dU}{dy};\\ y(\infty) = 0,\ y^{\prime}(0) = 0, \\U(y) = \frac{1}{4}y^{4}(\gamma +\alpha\ln^{2}y + \beta\ln^{4}y),$where ${}^{\prime}$ indicates differentiation with respect to $x$.The value of the constants in the equation are as follows: Mp = 2.435*10^18; \[Alpha] = 1.4*10^-5; \[Beta] = 6.3*10^-8; \[Gamma] = -0.013; \[Lambda]6 = 0; The appendix outlines the algorithm as follows, where Eq.(2.23) is referring to the differential equation (and boundary conditions) stated above. I am attempting to recreate this algorithm in Mathematica but my attempts thus far of simply trying numerically integrate the differential equation using a shooting method with boundary conditions specified in the Appendix has just yielded errors owing to finding complex infinities.Any help with this matter would be greatly appreciated. Attachments:
3 Replies
Sort By:
Posted 9 months ago
 Hello,I don't no if my code is correct!. Plot is smilar to your attachment. α = 1.4*10^-5; β = 6.3*10^-8; γ = -0.013; U = Re[D[1/4*y[x]^4*(γ + α*Log[y[x]]^2 + β*Log[y[x]]^4), y[x]]] sol = With[{ϵ = 1/10000}, NDSolve[{y''[x] + 3/x*y'[x] == U, y[5000] == 0, y'[ϵ] == 0}, y, {x, ϵ, 5000}, Method -> {"BoundaryValues" -> {"Shooting", "StartingInitialConditions" -> {y[ϵ] == 2, y'[ϵ] == 0}}}]] Plot[y[x] /. sol, {x, 1/10000, 200}, PlotRange -> All] Regards,MI
 Here is my updated code to see if this helps people. It still doesn't yield the required results but the algorithm follows the same procedure as in the paper. zero = 10^-10; inf = 10^9; ref = inf - 10^-3; \[Delta] = 0.01; DU = Re[D[ 1/4 y[x]^4*(\[Gamma] + \[Alpha]*(Log[y[x]])^2 + \[Beta]*(Log[ y[x]])^4), y[x]]]; eqn = y''[x] + 3 y'[x]/x == DU; B = 0.070; n = 0; k = 1; \[Epsilon] = 10^-10; \[CapitalDelta] = 1; \!$$\*OverscriptBox[\(\[CapitalDelta]$$, $$~$$]\) = 1; Bvalues = {}; absdiffvalues = {}; While[ \!$$\*OverscriptBox[\(\[CapitalDelta]$$, $$~$$]\) > \[Epsilon], eqn0 = B + B^3/8*(\[Gamma] + (\[Alpha]/2)* Log[B] + \[Alpha]*(Log[B])^2 + \[Beta]*(Log[ B])^3 + \[Beta]*(Log[B])^4) x^2; initialcon1 = eqn0 /. x -> zero; initialcon2 = D[eqn0, x] /. x -> zero; sol = NDSolve[{y''[x] + 3 y'[x]/x == DU, y[inf] == 0, y'[zero] == 0}, y[x], {x, zero, inf}, Method -> {"BoundaryValues" -> {"Shooting", "StartingInitialConditions" -> {y[zero] == initialcon1, y'[zero] == initialcon2}}}]; x2sol = x^2*y[x] /. (First@sol); \[CapitalDelta] = (x2sol /. x -> inf) - (x2sol /. x -> ref); \!$$\*OverscriptBox[\(\[CapitalDelta]$$, $$~$$]\) = Abs[\[CapitalDelta]]; AppendTo[absdiffvalues, \!$$\*OverscriptBox[\(\[CapitalDelta]$$, $$~$$]\)]; If[n >= 3, B0 = Part[Bvalues, n]; B1 = Part[Bvalues, n - 1]; B2 = Part[Bvalues, n - 2] ]; If[\[CapitalDelta] > \[Epsilon], If[B1 <= B0 <= B2 , k++; B += \[Delta]/(k*10); ];, B += \[Delta]/(k*10); ]; If[\[CapitalDelta] < -\[Epsilon], If[B2 <= B0 <= B1 , k++; B -= \[Delta]/(k*10); ];, B -= \[Delta]/(k*10); ]; AppendTo[Bvalues, B]; n++; If[Divisible[n, 500], Print[n]; Print[\[CapitalDelta]]; Print[B]; Print[k]; Print[initialcon1]; Print[initialcon2]; ]; ]; `Any help with this problem would still be greatly appreciated.