3
$\begingroup$

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,

  1. Is this a bug?
  2. Is there a work-around?
$\endgroup$
  • 1
    $\begingroup$ 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? $\endgroup$ – xzczd Jul 13 '17 at 2:53
  • 1
    $\begingroup$ @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. $\endgroup$ – bbgodfrey Jul 14 '17 at 0:06
  • $\begingroup$ @ilian You kindly answered my earlier question about parametric sensitivity. I would appreciate your thoughts about the behavior shown in the current question. Thanks. $\endgroup$ – bbgodfrey Jul 14 '17 at 1:16
  • 1
    $\begingroup$ 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. $\endgroup$ – ilian Jul 24 '17 at 14:58
  • 2
    $\begingroup$ Bug still present in MMA 11.2.0 $\endgroup$ – user58955 Sep 19 '17 at 15:06

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.