# Help in finding a good guesses for ode

This problem solved in 1 by @bbgodfrey. However, when I changed the parameters, the ode becomes very stiff and I did not know which best guess I need to use for y0. I am looking for :

 ParallelTable[Quiet[p[Z, g, k2, 1.5]], {g, 0.01, 0.03, .02}, {Z, 1000, 4000, 500},{k2, 0.0002, 100, 50}]


Does anyone know a method to estimate a good initial guess? I do not want to just put a random tries every time. I tried -.1 < y0 < -.08,but did not find a solution.

Here is my code:

p[Z0_, g0_, k0_, R0_] := Block[{Z = Z0,
g = Rationalize[g0, 0], k2 = Rationalize[k0, 0], ϵ = 10^-4, R = Rationalize[R0, 0]},
ps = ParametricNDSolveValue[{y''[r] + 2 y'[r]/r == k2 Sinh[y[r]], y[ϵ] == y0, y'[ϵ] == 0,
WhenEvent[r == 1, y'[r] -> y'[r] + Z g]}, {y, y'[R]}, {r, ϵ], R}, {y0},
Method -> "StiffnessSwitching", WorkingPrecision -> 20];
sol = FindRoot[Last[ps[y0]], {y0, -1}, Evaluated -> False][[1, 2]];
tot = 4 π sol NIntegrate[r^2 Exp[-First[ps[sol]][r]], {r, ϵ, R}];
L = Z/g*tot]
DistributeDefinitions[p];
tab = ParallelTable[Quiet[p[Z, g, 0.0002, 1.5]], {Z, 500, 800, 50}, {g, 0.02, 0.06, .005}];
ListContourPlot[tab, DataRange -> {{.02, .06}, {500, 800}}, PlotLegends -> Automatic,
FrameLabel -> {g, Z}, ImageSize -> Large, LabelStyle -> {Bold, Black, 15}]


Simply stated, a "good" initial guess for y0 in FindRoot is one for which ParametricNDSolve can integrate the ODE in the question from r = ϵ to r = R without errors. For the parameter ranges used in question 210296, -1 or any other negative number roughly in the range {-5, 0} worked well. The present question, however, introduces a much wider set of parameters, some of which accommodate only tiny ranges for the initial guesses. However, if the desired guess can be bracketed in y0, then the bracketed region can systematically be shrunk until a valid guess is obtained.

Specifically, fy in the code block below begins by bifurcating, yc = (yl + yu)/2, the bracketing region, and checks by means of tst = First[ps[yc]]["Domain"][[1, 2]] to determine whether ParametricNDSolve is able to integrate the ODE all the way to R for y0 = yc. If so, fy has found a good y00, which it Throws to p and terminates. Otherwise, it takes advantage of the fact that y'[tst] is negative when yc is more negative than the band of good guesses and is positive when yc is more positive than the band of good guesses in order to reduce the size of the bracketing region by a factor of two, after which the process is iterated until a good guess is found. Typically, less than ten iterations are required to find y00.

p[Z0_, g0_, k0_, R0_] := Block[{Z = Z0, g = Rationalize[g0, 0], k2 = Rationalize[k0, 0],
ϵ = 10^-4, R = Rationalize[R0, 0], yl = -8, yu = 0, ps, fy, y00, sol},
ps = ParametricNDSolveValue[{y''[r] + 2 y'[r]/r == k2 Sinh[y[r]], y[ϵ]] == y0,
y'[ϵ] == 0, WhenEvent[r == 1, y'[r] -> y'[r] + Z g]}, {y, y'}, {r, ϵ, R},
{y0}, Method -> "StiffnessSwitching", WorkingPrecision -> 20];
fy := (Do[yc = (yl + yu)/2;
tst = First[ps[yc]]["Domain"][[1, 2]];
If[tst == R, Throw[yc]];
If[Last[ps[yc]][tst] > 0, yu = yc, yl = yc], {i, 50}]; yc);
y00 = Which[R == First[ps[yl]]["Domain"][[1, 2]], yl,
R == First[ps[yu]]["Domain"][[1, 2]], yu, True, Catch[fy]];
sol = FindRoot[Last[ps[y0]][R], {y0, y00}, Evaluated -> False][[1, 2]];
L = Z/g 4 π sol NIntegrate[r^2 Exp[-First[ps[sol]][r]], {r, ϵ, R}]]


Applying p to the parameter set introduced by the present question yields.

DistributeDefinitions[p];
ParallelTable[Quiet[p[Z, g, k2, 1.5]], {Z, 1000, 4000, 500}, {g, 0.01, 0.03, .02},
{k2, 0.0002, 100, 50}]

(* {{{-1.36776*10^12, -14119.8}, {-1.422*10^12, -16956.}},
{{-3.13388*10^12, -33652.5}, {-3.2242*10^12, -40906.5}},
{{-5.62603*10^12, -63057.8}, {-5.79668*10^12, -75910.8}},
{{-8.84774*10^12, -102510.}, {-9.06488*10^12, -122341.}},
{{-1.2798*10^13, -152604.}, {-9.58419*10^10, -179866.}},
{{-1.74782*10^13, -213431.}, {-6.31954*10^10, -248535.}},
{{-2.2886*10^13, -287387.}, {-5.70976*10^10, -327845.}}} *)


The computation requires about ninety seconds on my six-processor computer, with only a small portion of that time used to obtain initial guesses. Note that the iteration over Z, which has most elements, is placed at the leftmost position to maximize the advantage of parallelization. Placing it as in the question increases the runtime to four minutes.

In closing, I would make two points. First, the "guess finder function", fy, fails if {yl, yu} does not bracket at least part of the band of good guesses. If fy does not converge to a good guess, sufficiently decreasing yl should solve the problem. Second, WorkingPrecision may need to be increased for some parameters to improve accuracy.

• Nice work @bbgodfrey, Why you need to use WhenEvent[Abs[y[r]] > 10^5, "StopIntegration"]? – aluuzz Dec 4 '19 at 21:28
• I really appreciate what have you done @bbgodfrey. I tried to understand what you have done here, but I did not figure out it. Could you please explain what is the name of iteration and why yl+yu/2. I am sorry for this question but I have no experience in iterative methods : fy := (Do[yc = (yl + yu)/2; tst = First[ps[yc]]["Domain"][[1, 2]]; If[tst == R, Throw[yc]]; If[Last[ps[yc]][tst] > 0, yu = yc, yl = yc], {i, 50}]; yc); y00 = Which[R == First[ps[yl]]["Domain"][[1, 2]], yl, R == First[ps[yu]]["Domain"][[1, 2]], yu, True, Catch[fy]]; – aluuzz Dec 5 '19 at 0:54
• @aluuzz No reason to be sorry for your question. I added an explanation to the answer. By the way, the second WhenEvent, which you asked about earlier, turned out not to have much value, so I deleted it. I also listed the variables used inside Block in its internal variables list, so that they would not interfere with any variables with the same names outside p. Typically, not necessary but good practice. – bbgodfrey Dec 5 '19 at 15:12