Message Boards Message Boards

0
|
10103 Views
|
5 Replies
|
0 Total Likes
View groups...
Share
Share this post:

[?] Speed up FindRoot with NIntegrate?

James Smith
Posted 8 years ago

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,
POSTED BY: James Smith
5 Replies
Sort By:
James Smith
James Smith
Posted 8 years ago

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.

POSTED BY: James Smith
Updating Name
Updating Name
Posted 8 years ago

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};
POSTED BY: Updating Name

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*)
POSTED BY: Mariusz Iwaniuk
James Smith
James Smith
Posted 8 years ago

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 http://support.wolfram.com/kb/12485 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.

POSTED BY: James Smith

First, please make sure you read and understand this:

http://support.wolfram.com/kb/12502

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:

http://support.wolfram.com/kb/12485

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

http://reference.wolfram.com/language/ref/ParametricNDSolveValue.html

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.
POSTED BY: Sean Clarke
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Tag limit exceeded
Note: Only the first five people you tag will receive an email notification; the other tagged names will appear as links to their profiles.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract