Skip to main content

numerics - How to implement custom NIntegrate integration strategies?


How can new integration strategies algorithms be used with NIntegrate?


This is a different type of extension than the extensions with new integration rules, as described in the answer for the question "How to implement custom integration rules for use by NIntegrate?". (Integration strategies perform and guide the core integration process, using or leveraging different integration rules and/or preprocessing algorithms.)



Answer



Motivation (for a new semi-symbolic integration strategy)



Consider the following integral, which cannot be done neigther by Integrate:


Integrate[BesselJ[y, x^3], {x, 0, ∞}, {y, 0, 1}]

(* Integrate[If[Re[y] > -(1/3), Gamma[1/6 + y/2]/(3*2^(2/3)*Gamma[5/6 + y/2]), 
     Integrate[BesselJ[y, x^3], {x, 0, Infinity}, 
   Assumptions -> Re[y] <= -(1/3)]], {y, 0, 1}] *)

nor NIntegrate:


NIntegrate[BesselJ[y, x^3], {x, 0, ∞}, {y, 0, 1}, 
 Method -> {"GlobalAdaptive", "MaxErrorIncreases" -> 2000}]



NIntegrate::slwcon: Numerical integration converging too slowly;


NIntegrate::eincr: The global error of the strategy GlobalAdaptive has ...



(* 0.524338 *)

Here is a plot of the integrand function over a much smaller domain:


Plot3D[BesselJ[y, x^3], {x, 0, 10}, {y, 0, 1}, PlotPoints -> {100, 10}, 
 MaxRecursion -> 5, PlotRange -> All, BoxRatios -> {10, 3}]



Because of the oscillatory nature of the integrand we can see why NIntegrate has difficulties.


On the other hand, Integrate can find the value of the integral integrating over $x$:


In[14]:= Integrate[BesselJ[y, x^3], {x, 0, ∞}]

Out[14]= ConditionalExpression[Gamma[1/6 + y/2]/(3*2^(2/3)*Gamma[5/6 + y/2]), 
 Re[y] > -(1/3)]

but it has problems integrating over $y$:



In[15]:= Integrate[BesselJ[y, x^3], {y, 0, 1}]

Out[15]= Integrate[BesselJ[y, x^3], {y, 0, 1}]

Since Integrate can do partially the integral along one of the axes, we can just take that symbolic expression and give it to NIntegrate for integration over the other axis.


Semi-symbolic NIntegrate implementation


Here we make an integration strategy that combines Integrate and NIntegrate -- it uses Integrate over some of the integration range(s) and then NIntegrate for the rest of the range(s) with the symbolic expressions obtained by Integrate.


The following defintion is for the initialization of the integration strategy SemiSymbolic.


Clear[SemiSymbolic];
Options[SemiSymbolic] = {"AnalyticalVariables" -> {}};

SemiSymbolicProperties = Options[SemiSymbolic][[All, 1]];
SemiSymbolic /: 
  NIntegrate`InitializeIntegrationStrategy[SemiSymbolic, nfs_, ranges_, 
   strOpts_, allOpts_] :=
  Module[{t, anVars},
   t = NIntegrate`GetMethodOptionValues[SemiSymbolic, SemiSymbolicProperties, 
     strOpts];
   If[t === $Failed, Return[$Failed]];
   {anVars} = t;
   SemiSymbolic[First /@ ranges, anVars]

   ];

This is the implementation of the integaration strategy SemiSymbolic:


SemiSymbolic[vars_, anVars_]["Algorithm"[regions_, opts___]] :=
  Module[{ranges, anRanges, funcs, t},

   ranges = Map[Flatten /@ Transpose[{vars, #@"OriginalBoundaries"}] &, regions];
   ranges = Map[Flatten, ranges, {-2}];
   anRanges = Map[Select[#, MemberQ[anVars, #[[1]]] &] &, ranges];
   ranges = Map[Select[#, ! MemberQ[anVars, #[[1]]] &] &, ranges];

   funcs = (#@"Integrand"[])@"FunctionExpression"[] & /@ regions;

   t = MapThread[
     Integrate[#1, Sequence @@ #2, 
       Assumptions -> (#[[2]] <= #[[1]] <= #[[3]] & /@ #3)] &, {funcs, 
      anRanges, ranges}];
   Print["SemiSymbolic::Integrate's result:", t];

   If[! FreeQ[t, Integrate], Return[$Failed]];


   Total[MapThread[
     NIntegrate[#1, Sequence @@ #2 // Evaluate, 
       Sequence @@ DeleteCases[opts, Method -> _] // Evaluate] &, {t, ranges}]]
   ];

(Note the implementation prints the intermediate result obtained by Integrate.)


Signatures


Initialization


We can see that the new rule SemiSymbolic is defined through TagSetDelayed for SemiSymbolic and NIntegrate`InitializeIntegrationStrategy. The rest of the arguments are:


nfs -- numerical function objects; several might be given depending on the integrand and ranges;



ranges -- a list of ranges for the integration variables;


strOpts -- the options given to the strategy;


allOpts -- all options given to NIntegrate.


Algorithm


StrategySymbol[strategyData___]["Algorithm"[regions_, opts___]] := ...

The algorithm can use regions objects as described in this answer of "Determining which rule NIntegrate selects automatically".


Remarks



Testing SemiSymbolic



The strategy works without (observable) problems for the motivational integral:


In[85]:= NIntegrate[BesselJ[y, x^3], {x, 0, Infinity}, {y, 0, 1}, 
 Method -> {SemiSymbolic, "AnalyticalVariables" -> {x}}]

During evaluation of In[85]:= SemiSymbolic::Integrate's result:{(2^-y HypergeometricPFQ[{1/6+y/2},{7/6+y/2,1+y},-(1/4)])/((1+3 y) Gamma[1+y])}

Out[85]= 0.524448

Note the printout for the intermediate result by Integrate.


Since SemiSymbolic passes inside its body the non-method NIntegrate options it was invoked with we can also see the sampling points used by SemiSymbolic through EvaluationMonitor.



res = 
 Reap@NIntegrate[BesselJ[y, x^3], {x, 0, Infinity}, {y, 0, 1}, 
   Method -> {SemiSymbolic, "AnalyticalVariables" -> {x}}, 
   EvaluationMonitor :> Sow[{x, y}]]

During evaluation of In[78]:= SemiSymbolic::Integrate's result:{(2^-y HypergeometricPFQ[{1/6+y/2},{7/6+y/2,1+y},-(1/4)])/((1+3 y) Gamma[1+y])}

(* {0.524448, {{{x, 0.00795732}, {x, 0.0469101}, {x, 0.122917}, {x, 
      0.230765}, {x, 0.360185}, {x, 0.5}, {x, 0.639815}, {x, 
      0.769235}, {x, 0.877083}, {x, 0.95309}, {x, 0.992043}, {x, 

      0.00397866}, {x, 0.023455}, {x, 0.0614583}, {x, 0.115383}, {x, 
      0.180092}, {x, 0.25}, {x, 0.319908}, {x, 0.384617}, {x, 
      0.438542}, {x, 0.476545}, {x, 0.496021}, {x, 0.503979}, {x, 
      0.523455}, {x, 0.561458}, {x, 0.615383}, {x, 0.680092}, {x, 
      0.75}, {x, 0.819908}, {x, 0.884617}, {x, 0.938542}, {x, 
      0.976545}, {x, 0.996021}}}} *)

ListPlot[res[[2, 1, All, 2]], Frame -> True]

enter image description here



Further tests


Below are some other tests / examples.


In[50]:= NIntegrate[x^2 + y^2 + z^2, {x, 0, 1}, {y, 0, 1}, {z, 0, 1}, 
 Method -> {SemiSymbolic, "AnalyticalVariables" -> {x, y}}]

During evaluation of In[50]:= SemiSymbolic::Integrate's result:{2/3+z^2}

Out[50]= 1.

Note that the symbolic integration was done over two variables.



Let us use the same integrand but with different range boundaries for the different variables in order to evaluate better the variable correspondence in the 2D sampling points pattern.


In[66]:= res = 
 Reap@NIntegrate[x^2 + y^2 + z^2, {x, 0, 1}, {y, 0, 2}, {z, 0, 10}, 
   Method -> {SemiSymbolic, "AnalyticalVariables" -> {x}}, 
   EvaluationMonitor :> Sow[{y, z}]]

During evaluation of In[66]:= SemiSymbolic::Integrate's result:{1/3+y^2+z^2}

Out[66]= {700., {{{1., 5.}, {1.35857, 5.}, {0.641431, 5.}, {1.94868, 
    5.}, {0.0513167, 5.}, {1., 6.79284}, {1., 3.20716}, {1., 

    9.74342}, {1., 0.256584}, {1.94868, 9.74342}, {1.94868, 
    0.256584}, {0.0513167, 0.256584}, {0.0513167, 9.74342}, {0.311753,
     1.55876}, {0.311753, 8.44124}, {1.68825, 1.55876}, {1.68825, 
    8.44124}}}}

In[69]:= ListPlot[res[[2, 1]], Frame -> True]

enter image description here


Another example


This Lebesgue integration implementation, AdaptiveNumericalLebesgueIntegration.m -- discussed in detail in "Adaptive numerical Lebesgue integration by set measure estimates" -- has implementations of integration strategy (and rules) with the complete signatures for the plug-in mechanism.



Comments

Post a Comment

Popular posts from this blog

plotting - Filling between two spheres in SphericalPlot3D

Manipulate[ SphericalPlot3D[{1, 2 - n}, {θ, 0, Pi}, {ϕ, 0, 1.5 Pi}, Mesh -> None, PlotPoints -> 15, PlotRange -> {-2.2, 2.2}], {n, 0, 1}] I cant' seem to be able to make a filling between two spheres. I've already tried the obvious Filling -> {1 -> {2}} but Mathematica doesn't seem to like that option. Is there any easy way around this or ... Answer There is no built-in filling in SphericalPlot3D . One option is to use ParametricPlot3D to draw the surfaces between the two shells: Manipulate[ Show[SphericalPlot3D[{1, 2 - n}, {θ, 0, Pi}, {ϕ, 0, 1.5 Pi}, PlotPoints -> 15, PlotRange -> {-2.2, 2.2}], ParametricPlot3D[{ r {Sin[t] Cos[1.5 Pi], Sin[t] Sin[1.5 Pi], Cos[t]}, r {Sin[t] Cos[0 Pi], Sin[t] Sin[0 Pi], Cos[t]}}, {r, 1, 2 - n}, {t, 0, Pi}, PlotStyle -> Yellow, Mesh -> {2, 15}]], {n, 0, 1}]
Post a Comment
Read more

plotting - Plot 4D data with color as 4th dimension

I have a list of 4D data (x position, y position, amplitude, wavelength). I want to plot x, y, and amplitude on a 3D plot and have the color of the points correspond to the wavelength. I have seen many examples using functions to define color but my wavelength cannot be expressed by an analytic function. Is there a simple way to do this? Answer Here a another possible way to visualize 4D data: data = Flatten[Table[{x, y, x^2 + y^2, Sin[x - y]}, {x, -Pi, Pi,Pi/10}, {y,-Pi,Pi, Pi/10}], 1]; You can use the function Point along with VertexColors . Now the points are places using the first three elements and the color is determined by the fourth. In this case I used Hue, but you can use whatever you prefer. Graphics3D[ Point[data[[All, 1 ;; 3]], VertexColors -> Hue /@ data[[All, 4]]], Axes -> True, BoxRatios -> {1, 1, 1/GoldenRatio}]
Post a Comment
Read more

plotting - Adding a thick curve to a regionplot

Suppose we have the following simple RegionPlot: f[x_] := 1 - x^2 g[x_] := 1 - 0.5 x^2 RegionPlot[{y < f[x], f[x] < y < g[x], y > g[x]}, {x, 0, 2}, {y, 0, 2}] Now I'm trying to change the curve defined by $y=g[x]$ into a thick black curve, while leaving all other boundaries in the plot unchanged. I've tried adding the region $y=g[x]$ and playing with the plotstyle, which didn't work, and I've tried BoundaryStyle, which changed all the boundaries in the plot. Now I'm kinda out of ideas... Any help would be appreciated! Answer With f[x_] := 1 - x^2 g[x_] := 1 - 0.5 x^2 You can use Epilog to add the thick line: RegionPlot[{y < f[x], f[x] < y < g[x], y > g[x]}, {x, 0, 2}, {y, 0, 2}, PlotPoints -> 50, Epilog -> (Plot[g[x], {x, 0, 2}, PlotStyle -> {Black, Thick}][[1]]), PlotStyle -> {Directive[Yellow, Opacity[0.4]], Directive[Pink, Opacity[0.4]],
Post a Comment
Read more