# Unexpected Behavior of Parametric Sensitivity in ParametricNDSolveValue

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.000115.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,

1. Is this a bug?
2. Is there a work-around?
• In v9.0.1, "Domain" seems not to change in the first case: i.stack.imgur.com/EqKkt.png and only changes slightly in the second case: 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 2:53
• @xzczd Thank you for your advice. I tried 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. – bbgodfrey Jul 14 '17 at 0:06
• @ilian You kindly answered my earlier question about parametric sensitivity. I would appreciate your thoughts about the behavior shown in the current question. Thanks. – bbgodfrey Jul 14 '17 at 1:16
• Thank you, now I can confirm the same behavior. A bug report has been opened, but I am not aware of a workaround for the moment. – ilian Jul 24 '17 at 14:58
• Bug still present in MMA 11.2.0 – user58955 Sep 19 '17 at 15:06