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: i.sstatic.net/EqKkt.png and only changes slightly in the second case: i.sstatic.net/1tPAO.png , I think it's a backslide at least. This might be related to this issue so try the workaround there? $\endgroup$Method -> {"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. $\endgroup$