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

mathematical optimization - Minimizing using indices, error: Part::pkspec1: The expression cannot be used as a part specification

I want to use Minimize where the variables to minimize are indices pointing into an array. Here a MWE that hopefully shows what my problem is. vars = u@# & /@ Range[3]; cons = Flatten@ { Table[(u[j] != #) & /@ vars[[j + 1 ;; -1]], {j, 1, 3 - 1}], 1 vec1 = {1, 2, 3}; vec2 = {1, 2, 3}; Minimize[{Total@((vec1[[#]] - vec2[[u[#]]])^2 & /@ Range[1, 3]), cons}, vars, Integers] The error I get: Part::pkspec1: The expression u[1] cannot be used as a part specification. >> Answer Ok, it seems that one can get around Mathematica trying to evaluate vec2[[u[1]]] too early by using the function Indexed[vec2,u[1]] . The working MWE would then look like the following: vars = u@# & /@ Range[3]; cons = Flatten@{ Table[(u[j] != #) & /@ vars[[j + 1 ;; -1]], {j, 1, 3 - 1}], 1 vec1 = {1, 2, 3}; vec2 = {1, 2, 3}; NMinimize[ {Total@((vec1[[#]] - Indexed[vec2, u[#]])^2 & /@ R...
Post a Comment
Read more

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...
Post a Comment
Read more

What is and isn't a valid variable specification for Manipulate?

I have an expression whose terms have arguments (representing subscripts), like this: myExpr = A[0] + V[1,T] I would like to put it inside a Manipulate to see its value as I move around the parameters. (The goal is eventually to plot it wrt one of the variables inside.) However, Mathematica complains when I set V[1,T] as a manipulated variable: Manipulate[Evaluate[myExpr], {A[0], 0, 1}, {V[1, T], 0, 1}] (*Manipulate::vsform: Manipulate argument {V[1,T],0,1} does not have the correct form for a variable specification. >> *) As a workaround, if I get rid of the symbol T inside the argument, it works fine: Manipulate[ Evaluate[myExpr /. T -> 15], {A[0], 0, 1}, {V[1, 15], 0, 1}] Why this behavior? Can anyone point me to the documentation that says what counts as a valid variable? And is there a way to get Manpiulate to accept an expression with a symbolic argument as a variable? Investigations I've done so far: I tried using variableQ from this answer , but it says V[1...
Post a Comment
Read more