Bug introduced in 10.4 or earlier and continuing through 11.3
Submitted as CASE:3916971
While exploring alternative methods of solving 33538, I encountered difficulties with the parametric sensitivity capability of ParametricNDSolveValue. (Useful background on parametric sensitivity is provided here and here.) Specifically, for
r0 = 1/10000; inf = 50; ξ = 15/10; n = 1;
pa = 172416/390625; pf = 5298661942732265830186933/7730236542592271015625000;
numSol = ParametricNDSolveValue[
{r D[1/r D[A[r], r], r] - ξ^2 F[r]^2 (A[r] - 1) == 0,
1/r D[r D[F[r], r], r] - n^2/r^2 F [r] (A[r] - 1)^2 - 1/2 F[r] (F[r]^2 - 1) == 0,
A[r0] == a r0^2 - If[n == 1, f^2 r0^n ξ^2/8, 0],
F[r0] == f r0^n (1 - (1 + 4 a n^2) r0^2/(8 (n + 1))),
A'[r0] == 2 a r0 - If[n == 1, r0^3 f^2 ξ^2/2, 0],
F'[r0] == f r0^(n - 1) (n - (2 + n) (1 + 4 a n^2) r0^2/(8 (1 + n))),
WhenEvent[A'[r] < 0 || F'[r] < 0 || A[r] > 1 || F[r] > 1, "StopIntegration"]},
{A, F}, {r, r0, inf}, {a, f}];
the domain of integration changes when the sensitivity of the solution to infinitesimal variation of the values of the parameters is determined.
numSol[a, f] /. {a -> pa, f -> pf};
(First@%)["Domain"][[1, 2]]
D[numSol[a, f], a] /. {a -> pa, f -> pf};
(First@%)["Domain"][[1, 2]]
D[numSol[a, f], f] /. {a -> pa, f -> pf};
(First@%)["Domain"][[1, 2]]
numSol[a, f] /. {a -> pa, f -> pf};
(First@%)["Domain"][[1, 2]]
(* 6.21431 *)
(* 6.05906 *)
(* 6.05906 *)
(* 6.05906 *)
This change does not occur, if the sensitivity determination is not performed. To see this, again execute the code defining numSol, and then execute
numSol[a, f] /. {a -> pa, f -> pf}
(First@%)["Domain"][[1, 2]]
numSol[a, f] /. {a -> pa, f -> pf}
(First@%)["Domain"][[1, 2]]
(* 6.21431 *)
(* 6.21431 *)
This strange behavior is exacerbated, when WorkingPrecision is set explicitly. For instance,
numSol = ParametricNDSolveValue[
{r D[1/r D[A[r], r], r] - ξ^2 F[r]^2 (A[r] - 1) == 0,
1/r D[r D[F[r], r], r] - n^2/r^2 F [r] (A[r] - 1)^2 - 1/2 F[r] (F[r]^2 - 1) == 0,
A[r0] == a r0^2 - If[n == 1, f^2 r0^n ξ^2/8, 0],
F[r0] == f r0^n (1 - (1 + 4 a n^2) r0^2/(8 (n + 1))),
A'[r0] == 2 a r0 - If[n == 1, r0^3 f^2 ξ^2/2, 0],
F'[r0] == f r0^(n - 1) (n - (2 + n) (1 + 4 a n^2) r0^2/(8 (1 + n))),
WhenEvent[A'[r] < 0 || F'[r] < 0 || A[r] > 1 || F[r] > 1, "StopIntegration"]},
{A, F}, {r, r0, inf}, {a, f}, WorkingPrecision -> $MachinePrecision];
in which case the calls to numSol listed above return
(* 6.214311927593723 *)
(* 7.462342033621696 *)
ParametricNDSolveValue::ndsz: At r == 0.0001`15.954589770191003, step size is effectively zero; singularity or stiff system suspected.
(* 0.0001000000000000000 *)
(* 0.0001000000000000000 *)
The computations were performed with version 11.1.1 and also with 10.4.1 on Windows 10 (64 bit). My specific questions are,
- Is this a bug?
- Is there a work-around?
"Domain"seems not to change in the first case: https://i.stack.imgur.com/EqKkt.png and only changes slightly in the second case: https://i.stack.imgur.com/1tPAO.png , I think it's a backslide at least. This might be related to this issue so try the workaround there? – xzczd Jul 13 '17 at 02:53Method -> {"DiscontinuityProcessing" -> False}, as you suggested, but it did not help. Also, the bug described in the linked question purportedly has been fixed in 11.1.1. I plan to submit this as a bug. – bbgodfrey Jul 14 '17 at 00:06pf, now corrected. (I had dropped the first digit in the numerator.) Please try again to reproduce my result. Probably, the value ofpfmatters a lot, because this is a separatrix computation. The correct value ofpfallows integration far along the separatrix, where tiny changes in the initial conditions (and presumably anything else) can result in big changes to the output. In contrast, the earlier incorrect value ofpfscarcely allowed the integration to reach the separatrix, so any results it returned were fairly insensitive to initial conditions. – bbgodfrey Jul 19 '17 at 20:24