Skip to main content

numerics - How to implement custom integration rules for use by NIntegrate?


How can NIntegrate be extended with custom implementation of integration rules?


This answer of the question "Monte Carlo integration with random numbers generated from a Gaussian distribution" shows customization or a general rule ("MonteCarloRule"). This question is about defining new rules.



A related question is "How to implement custom NIntegrate integration strategies?".



Answer



The simplest way to make new NIntegrate algorithms is by user defined integration rules. Below are given examples using a simple rule (the Simpson rule) and how NIntegrate's framework can utilize the new rule implementations with its algorithms. (Adaptive, symbolic processing, and singularity handling algorithms are seamlessly applied.)


Basic 1D rule implementation (Simpson rule)


This easiest way to add a new rule is to use NIntegrate's symbol GeneralRule. Such a rule is simply initialized with a list of three elements:


{abscissas, integral weights, error weights}

The Simpson rule:


$$ \int_0^1 f(x)dx \approx \frac{1}{6} \lgroup f(0)+4 f(\frac{1}{2})+f(1) \rgroup$$


is implemented with the following definition:



SimpsonRule /: 
 NIntegrate`InitializeIntegrationRule[SimpsonRule, nfs_, ranges_, ruleOpts_, allOpts_] := 
 NIntegrate`GeneralRule[{{0, 1/2, 1}, {1/6, 4/6, 1/6}, {1/6, 4/6, 1/6} - {1/2, 0, 1/2}}]

The error weights are calculated as the difference between the Simson rule and the trapezoidal rule.


Signature


We can see that the new rule SimpsonRule is defined through TagSetDelayed for SimpsonRule and NIntegrate`InitializeIntegrationRule. 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;


ruleOpts -- the options given to the rule;



allOpts -- all options given to NIntegrate.


Note that here we discuss the rule algorithm initialization only. The discussed intializations produce general rules for which there is an implemented computation algorithm.


(Explaining the making of definitions for integration rule computation algorithms is postponed for now. See the MichaelE2 answer or this blog post and related package AdaptiveNumericalLebesgueIntegration.m for examples of how to hook-up integration rules computation algorithms.)


NIntegrate's plug-in mechanism is fairly analogous to NDSolve's -- see the tutorial "NDSolve Method Plugin Framework".


Basic 1D rule tests


Here is the test with the SimpsonRule implemented above.


NIntegrate[Sqrt[x], {x, 0, 1}, Method -> SimpsonRule]

(* 0.666667 *)


Here are the sampling points of the integration above.


k = 0;
ListPlot[Reap[
   NIntegrate[Sqrt[x], {x, 0, 1}, Method -> SimpsonRule, 
    EvaluationMonitor :> Sow[{x, ++k}]]][[2, 1]], 
 PlotTheme -> "Detailed", ImageSize -> Large]

enter image description here


Multi-panel Simpson rule implementation


Here is an implementation of the multi-panel Simson rule:



Options[MultiPanelSimpsonRule] = {"Panels" -> 5};
MultiPanelSimpsonRuleProperties = Part[Options[MultiPanelSimpsonRule], All, 1];
MultiPanelSimpsonRule /: 
  NIntegrate`InitializeIntegrationRule[MultiPanelSimpsonRule, nfs_, ranges_, 
   ruleOpts_, allOpts_] :=

  Module[{t, panels, pos, absc, weights, errweights},

   t = NIntegrate`GetMethodOptionValues[MultiPanelSimpsonRule, 
     MultiPanelSimpsonRuleProperties, ruleOpts];

   If[t === $Failed, Return[$Failed]];
   {panels} = t;

   If[! TrueQ[NumberQ[panels] && 1 <= panels < Infinity], 
    pos = NIntegrate`OptionNamePosition[ruleOpts, "Panels"];
    Message[NIntegrate::intpm, ruleOpts, {pos, 2}];
    Return[$Failed];
    ];

   weights = Table[{1/6, 4/6, 1/6}, {panels}];

   weights = 
    Fold[Join[Drop[#1, -1], {#1[[-1]] + #2[[1]]}, Rest[#2]] &, First[weights],
       Rest[weights]]/panels;
   {absc, errweights, t} = 
    NIntegrate`TrapezoidalRuleData[(Length[weights] + 1)/2, 
     WorkingPrecision /. allOpts];
   NIntegrate`GeneralRule[{absc, weights, (weights - errweights)}]
   ];

Multi-panel Simpson rule tests



Here is an integral calculation with the multi-panel Simson rule


NIntegrate[Sqrt[x], {x, 0, 1}, 
 Method -> {MultiPanelSimpsonRule, "Panels" -> 12}]

(* 0.666667 *)

Here are the sampling points of the integration above:


k = 0;
ListPlot[Reap[
   NIntegrate[Sqrt[x], {x, 0, 1}, 

    Method -> {MultiPanelSimpsonRule, "Panels" -> 12}, MaxRecursion -> 10, 
    EvaluationMonitor :> Sow[{x, ++k}]]][[2, 1]]]

enter image description here


Note the traces of the "DoubleExponential" singularity handler application on the right side around 220th and 750th sampling points.


Two dimensional integration with a Cartisian product of Simpson multi-panel rules


The 1D multi-panel rule implemented above can be used for multi-dimensional integration.


This is what we get with NIntegrate's default method


NIntegrate[Sqrt[x + y], {x, 0, 1}, {y, 0, 1}]


(* 0.975161 *)

Here is the estimate with the custom multi-panel rule:


NIntegrate[Sqrt[x + y], {x, 0, 1}, {y, 0, 1}, 
 Method -> {MultiPanelSimpsonRule, "Panels" -> 5}, MaxRecursion -> 10]

(* 0.975161 *)

Note that the command above is equivalent to:


NIntegrate[Sqrt[x + y], {x, 0, 1}, {y, 0, 1}, 

 Method -> {"CartesianRule", 
   Method -> {MultiPanelSimpsonRule, "Panels" -> 5}}, MaxRecursion -> 10]

(* 0.975161 *)

Here is a plot of the sampling points:


k = 0;
ListPlot[Reap[
   NIntegrate[Sqrt[x + y], {x, 0, 1}, {y, 0, 1}, PrecisionGoal -> 5, 
    Method -> {MultiPanelSimpsonRule, "Panels" -> 5}, MaxRecursion -> 10, 

    EvaluationMonitor :> Sow[{x, y}]]][[2, 1]]]

enter image description here


Note the trace of the application of the singularity handler "DuffyCoordinates" at the left-bottom corner.


A complete example


This Lebesgue integration implementation, AdaptiveNumericalLebesgueIntegration.m -- discussed in detail in "Adaptive numerical Lebesgue integration by set measure estimates" -- has implementations of integration 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 - Mathematica: 3D plot based on combined 2D graphs

I have several sigmoidal fits to 3 different datasets, with mean fit predictions plus the 95% confidence limits (not symmetrical around the mean) and the actual data. I would now like to show these different 2D plots projected in 3D as in but then using proper perspective. In the link here they give some solutions to combine the plots using isometric perspective, but I would like to use proper 3 point perspective. Any thoughts? Also any way to show the mean points per time point for each series plus or minus the standard error on the mean would be cool too, either using points+vertical bars, or using spheres plus tubes. Below are some test data and the fit function I am using. Note that I am working on a logit(proportion) scale and that the final vertical scale is Log10(percentage). (* some test data *) data = Table[Null, {i, 4}]; data[[1]] = {{1, -5.8}, {2, -5.4}, {3, -0.8}, {4, -0.2}, {5, 4.6}, {1, -6.4}, {2, -5.6}, {3, -0.7}, {4, 0.04}, {5, 1.0}, {1, -6.8}, {2, -4.7}, {3, -1.
Post a Comment
Read more