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

Popular posts from this blog

front end - keyboard shortcut to invoke Insert new matrix

I frequently need to type in some matrices, and the menu command Insert > Table/Matrix > New... allows matrices with lines drawn between columns and rows, which is very helpful. I would like to make a keyboard shortcut for it, but cannot find the relevant frontend token command (4209405) for it. Since the FullForm[] and InputForm[] of matrices with lines drawn between rows and columns is the same as those without lines, it's hard to do this via 3rd party system-wide text expanders (e.g. autohotkey or atext on mac). How does one assign a keyboard shortcut for the menu item Insert > Table/Matrix > New... , preferably using only mathematica? Thanks! Answer In the MenuSetup.tr (for linux located in the $InstallationDirectory/SystemFiles/FrontEnd/TextResources/X/ directory), I changed the line MenuItem["&New...", "CreateGridBoxDialog"] to read MenuItem["&New...", "CreateGridBoxDialog", MenuKey["m", Modifiers-...

How to thread a list

I have data in format data = {{a1, a2}, {b1, b2}, {c1, c2}, {d1, d2}} Tableform: I want to thread it to : tdata = {{{a1, b1}, {a2, b2}}, {{a1, c1}, {a2, c2}}, {{a1, d1}, {a2, d2}}} Tableform: And I would like to do better then pseudofunction[n_] := Transpose[{data2[[1]], data2[[n]]}]; SetAttributes[pseudofunction, Listable]; Range[2, 4] // pseudofunction Here is my benchmark data, where data3 is normal sample of real data. data3 = Drop[ExcelWorkBook[[Column1 ;; Column4]], None, 1]; data2 = {a #, b #, c #, d #} & /@ Range[1, 10^5]; data = RandomReal[{0, 1}, {10^6, 4}]; Here is my benchmark code kptnw[list_] := Transpose[{Table[First@#, {Length@# - 1}], Rest@#}, {3, 1, 2}] &@list kptnw2[list_] := Transpose[{ConstantArray[First@#, Length@# - 1], Rest@#}, {3, 1, 2}] &@list OleksandrR[list_] := Flatten[Outer[List, List@First[list], Rest[list], 1], {{2}, {1, 4}}] paradox2[list_] := Partition[Riffle[list[[1]], #], 2] & /@ Drop[list, 1] RM[list_] := FoldList[Transpose[{First@li...

functions - Get leading series expansion term?

Given a function f[x] , I would like to have a function leadingSeries that returns just the leading term in the series around x=0 . For example: leadingSeries[(1/x + 2)/(4 + 1/x^2 + x)] x and leadingSeries[(1/x + 2 + (1 - 1/x^3)/4)/(4 + x)] -(1/(16 x^3)) Is there such a function in Mathematica? Or maybe one can implement it efficiently? EDIT I finally went with the following implementation, based on Carl Woll 's answer: lds[ex_,x_]:=( (ex/.x->(x+O[x]^2))/.SeriesData[U_,Z_,L_List,Mi_,Ma_,De_]:>SeriesData[U,Z,{L[[1]]},Mi,Mi+1,De]//Quiet//Normal) The advantage is, that this one also properly works with functions whose leading term is a constant: lds[Exp[x],x] 1 Answer Update 1 Updated to eliminate SeriesData and to not return additional terms Perhaps you could use: leadingSeries[expr_, x_] := Normal[expr /. x->(x+O[x]^2) /. a_List :> Take[a, 1]] Then for your examples: leadingSeries[(1/x + 2)/(4 + 1/x^2 + x), x] leadingSeries[Exp[x], x] leadingSeries[(1/x + 2 + (1 - 1/x...