[?] Speed up FindRoot with NIntegrate?

Questions
,
1

Dear all,
I am trying to find the root of a function F which is defined through a numerical integral with NIntegrate:

F[\[Mu]C_] := 
 1/NN NIntegrate[
    r Exp[-((NN - 1 )/2) r^2 + 
       NN Log[ BesselJ[0, Sqrt[-2 \[Mu]C] r]]] ((
        NN BesselJ[1, Sqrt[-2 \[Mu]C] r] r)/(
        BesselJ[0, Sqrt[-2 \[Mu]C] r] Sqrt[-2 \[Mu]C])
         BesselJ[
         0, ((NN - 1) r)/Sqrt[-2 \[Mu]C] Sqrt[\[Mu]M.\[Mu]M]] - 
       BesselJ[1, ((NN - 1) r)/Sqrt[-2 \[Mu]C]
           Sqrt[\[Mu]M.\[Mu]M]] ((NN - 1) r)/(-2 \[Mu]C)^(3/2)
         Sqrt[\[Mu]M.\[Mu]M]), {r, 0, \[Infinity]}, 
    Method -> {"ExtrapolatingOscillatory"}]/
   NIntegrate[
    r Exp[-((NN - 1 )/2) r^2 + 
       NN Log[ BesselJ[0, Sqrt[-2 \[Mu]C] r]]] BesselJ[
      0, ((NN - 1) r)/Sqrt[-2 \[Mu]C] Sqrt[\[Mu]M.\[Mu]M]], {r, 
     0, \[Infinity]}, Method -> {"ExtrapolatingOscillatory"}]

I then set the parameters that appear in F

NN = 60;
\[Mu]M = {1, -1};`

and try to find a root of the following equation

FindRoot[F[\[Mu]] == -0.4, {\[Mu], -1, -4, 0}, Evaluated -> False] .

The root exists, but FindRoot keeps running for a long time > 10 mins without any output. Do you know how to fix this problem?

Thank you
Best,

2

First, please make sure you read and understand this:

The basic point I’d like to make is that your function, “F” should be made to only evaluate when given a numerical value.

Here are some basic suggestions on how to speed up the evaluation of NIntegrate:

Wolfram Support Quick Answers: How do I accelerate NIntegrate evaluations?

Another suggestion is to rewrite your Integrals as differential equations and then to use ParametricNDSolve. This might work better.

Lastly, there is no substitute for trying to break your Integrals down and find decent simplifications or at worse approximations. Approximations can be used to find a good approximate root at least. Your integrals have two important properties (1) They seem to go to 0 very quickly and (2) they are highly oscillatory.

(1) In this case, maybe you want to do a change of variable so that the integral is going from 0 to 1 instead of 0 to Infinity.
(2) You can sometimes “tame” highly oscillatory integrals with integration by parts. As an example, consider integrating Sin[1/x] E^x. The Sin[1/x] oscillates around 0 excessively. But the good news is that it has a symbolic (or at least relatively easy to evaluate) integral. Let’s call the integral of Sin[1/x] “intSinoverx”. Because integration is an inherently smoothing operation, “inSinoverx” is far less oscillatory than Sin[1/x] and we can rewrite our integral using it. You can find what this integral actually is using Mathematica. Using integration by parts we can rewrite the integral as:

intSinoverx*E^x - Integrate[inSinoverx*E^x,x]

Anyway, there’s an infinite number of things you can try.

3

First, I know the option _?NumericQ , and I did not include it intentionally in the minimal working example because it does not solve the problem, and complicates the code.

The problem here is rather that the function F can be computed relatively quickly when called individually with a given argument, while it takes much longer to be evaluated when called from FindRoot with the exact same argument, see this minimal working example , where I make the function print the argument with which it was called. On my laptop, the time lapse for the first call is ˜ 2 seconds, while the time lapse between the first and second call of the function in FindRoot is much larger.

Unfortunately the suggestions reported in Wolfram Support Quick Answers: How do I accelerate NIntegrate evaluations? do not solve this problem. Trying to manipulate the integral analytically could help, but one needs to understand the basic problem above before trying other things.

4

Let’s go to details:

NN = 60;
?M = {1, -1};
f[?C_?NumericQ] := 1/NN NIntegrate[r Exp[-((NN - 1)/2) r^2 + NN Log[BesselJ[0, Sqrt[-2 ?C] r]]] ((NN BesselJ[1, 
Sqrt[-2 ?C] r] r)/(BesselJ[0, Sqrt[-2 ?C] r] Sqrt[-2 ?C]) BesselJ[0, ((NN - 1) r)/Sqrt[-2 ?C] Sqrt[?M.?M]] - BesselJ[1, ((NN - 1) r)/
Sqrt[-2 ?C] Sqrt[?M.?M]] ((NN - 1) r)/(-2 ?C)^(3/2) Sqrt[?M.?M]), {r, 0, ?}, 
Method -> {"LevinRule", "LevinFunctions" -> {BesselJ}, "MaxOrder" -> 100}]/
NIntegrate[r Exp[-((NN - 1)/2) r^2 + NN Log[BesselJ[0, Sqrt[-2 ?C] r]]] BesselJ[0, ((NN - 1) r)/Sqrt[-2 ?C] Sqrt[?M.?M]], {r, 0, ?}, 
Method -> {"LevinRule", "LevinFunctions" -> {BesselJ}, "MaxOrder" -> 100}]
 
sol = FindRoot[f[?] + 4/10 == 0, {?, -1, -2}] // Quiet
(* {? -> -0.962644}*)
f[?] /. sol // Quiet
(*-0.4*)
5

You could use the Levin rule, with a precomputed symbolic analysis. NIntegrate`LevinIntegrandReduce is discussed in NIntegrate Integration Rules. For some reason li["Rules"] does not return rules that are all suitable for input to the "LevinRule" method, so they need some processing.

NN = 60;
?M = {1, -1};

ClearAll[int1, int2, li1, li2, dli1, dli2, liMeth, FLR, dFLR, F1, F2, dF1, dF2];
int1[?C_] := r Exp[-((NN - 1)/2) r^2 + NN Log[BesselJ[0, Sqrt[-2 ?C] r]]]
  ((NN BesselJ[1, Sqrt[-2 ?C] r] r)/(BesselJ[0, Sqrt[-2 ?C] r] Sqrt[-2 ?C]) *
      BesselJ[0, ((NN - 1) r)/Sqrt[-2 ?C] Sqrt[?M.?M]] - 
     BesselJ[1, ((NN - 1) r)/Sqrt[-2 ?C] Sqrt[?M.?M]] ((NN - 1) r)/(-2 ?C)^(3/2) Sqrt[?M.?M]);
int2[?C_] := 
  r Exp[-((NN - 1)/2) r^2 + NN Log[BesselJ[0, Sqrt[-2 ?C] r]]] BesselJ[0, ((NN - 1) r)/Sqrt[-2 ?C] Sqrt[?M.?M]];
li1[?C_] = NIntegrate`LevinIntegrandReduce[int1[?C], {r}];
li2[?C_] = NIntegrate`LevinIntegrandReduce[int2[?C], {r}];
liMeth[li_] := Flatten@{"LevinRule",
    Rest@li["Rules"] /. HoldPattern["DifferentialMatrices" -> dmm_] :>  "DifferentialMatrix" -> First@dmm,
    "MaxOrder" -> 100};
FLR[?C_?NumericQ] := 1/NN NIntegrate[int1[?C], {r, 0, ?}, Method -> liMeth[li1[?C]]]/
    NIntegrate[int2[?C], {r, 0, ?}, Method -> liMeth[li2[?C]]];

AbsoluteTiming[
 sol = FindRoot[FLR[m] == -0.4, {m, -1, -2}]
 ]
(*  {4.28894, {m -> -0.962643}}  *)

Check:

FLR[m] /. sol
(*  -0.4  *)

An alternative definition of liMeth could be

liMeth[li_] := {"LevinRule",
   "Kernel" -> li["Kernel"],
   "Amplitude" -> li@"Amplitude",
   "AdditiveTerm" -> li@"AdditiveTerm",
   "DifferentialMatrix" -> First@li@"DifferentialMatrices",
   "MaxOrder" -> 100};
6

Thank you for the very useful suggestions about decomposing identifying symbolically the oscillatory part of the integral with LevinRule and its options. I am still getting a few error messages, or warnings, with it, but it is definitely a good idea to try.

Also, I realized how to fix the problem that I had the minimal working example above, i.e., that function evaluations were fast when calling the function individually, but slow when calling the function from FindRoot: this seems to be due to an attempt to process the integral symbolically by FindRoot, and can be solved by inserting the option Method → {Automatic, “SymbolicProcessing” → 0} in NIntegrate.