Stopping NDSolve when encountering stiffness
NDSolve
stops itself when it believes that it has encountered stiffness. To take advantage of this, rewrite test
as
test = ParametricNDSolveValue[{(Derivative[2][t][r] + (2/r)*
Derivative[1][t][r] - D[Potential[x], x] /. x -> t[r]) == 0,
t[10^(-12)] == d, Derivative[1][t][10^(-12)] == 0}, t, {r, 10^(-12), 1}, {d}]
In other words, obtain the entire solution t
, not just t[.5]
. Then use, for instance,
test[.001]["Domain"][[1, 2]]
(* 0.000383406 *)
which shows that NDSolve
detected stiffness at x == 0.000383406
. On the other hand,
test[.002]["Domain"][[1, 2]]
(* 1. *)
shows that NDSolve
did not encounter stiffness throughout the range of integration. Applying this to the code in the question that seeks the value of d
for which stiffness first occurs gives
c = 0.0005; d = 0.002; Under = 1; Monitor[
Quiet[While[c >= 10^(-18), While[Under == 1, d = d - c; If[test[d]["Domain"][[1, 2]] < 1,
Under = 0];]; d = d + c; c = c/10; s = d; Under = 1;];], d]
NumberForm[s, 16]
(* 0.001701281449747991 *)
Note that c
is used here instead of C
, which is a reserved symbol in Mathematica.
Addendum: Faster Solution
The solution above, patterned after that in the question, requires 108 sec
on my PC to obtain an answer. The following, which brackets the solution and successively reduces the bounds by a factor of two on each iteration,
dl = c; du = d;
Do[dt = (dl + du)/2; If[Quiet@test[dt]["Domain"][[1, 2]] < 1, dl = dt, du = dt], {i, 20}]
NumberForm[du, 16]
(* 0.001701281449747996 *)
takes 0.000112592 sec
.