Function Repository Resource:

PaduaIntegrate

Source Notebook

Numerically integrate a function using the Padua points

Contributed by: Jan Mangaldan

ResourceFunction["PaduaIntegrate"][f,{x,xmin,xmax},{y,ymin,ymax}]

gives a numerical approximation to the multiple integral by evaluating f at the Padua points.

Details and Options

The Padua points are a unisolvent point set that can be used to interpolate and integrate a bivariate function over a rectangular domain. ResourceFunction["PaduaIntegrate"] evaluates f at the Padua points, and integrates the interpolating polynomial passing through them.
ResourceFunction["PaduaIntegrate"] takes the following options:
InterpolationOrder15order of the interpolating polynomial generated
"PaduaType"1type of Padua points to use
WorkingPrecisionMachinePrecisionthe precision used in internal computations
Possible types of Padua points are 1, 2, 3 and 4.

Examples

Basic Examples (2) 

Numerically integrate the function Exp[-x2-y2]:

In[1]:=
ResourceFunction["PaduaIntegrate"][
 Exp[-x^2 - y^2], {x, -1, 1}, {y, -1, 1}]
Out[1]=

Compare with the exact answer:

In[2]:=
% - \[Pi] Erf[1]^2
Out[2]=

Scope (3) 

A test function due to Franke:

In[3]:=
franke[x_, y_] := 3/4 Exp[-(((9 x - 2)^2 + (9 y - 2)^2)/4)] + 3/4 Exp[-((9 x + 1)^2/49) - (9 y + 1)/10] + 1/2 Exp[-(((9 x - 7)^2 + (9 y - 3)^2)/4)] - 1/5 Exp[-(9 x - 4)^2 - (9 y - 7)^2]

Integrate the function over a rectangular domain:

In[4]:=
ResourceFunction["PaduaIntegrate"][
 franke[x, y], {x, 0, 1}, {y, 0, 1/2}]
Out[4]=

Compare with the result of NIntegrate:

In[5]:=
% - NIntegrate[franke[x, y], {x, 0, 1}, {y, 0, 1/2}]
Out[5]=

Options (3) 

InterpolationOrder (1) 

Use a degree 25 interpolant for integrating the Dixon–Szegö function:

In[6]:=
ResourceFunction[
 "PaduaIntegrate"][(4 - 2.1 x^2 + x^4/3) x^2 + x y + 4 (y^2 - 1) y^2, {x, -2, 2}, {y, -1.25, 1.25}, InterpolationOrder -> 25]
Out[6]=

PaduaType (1) 

Use type-3 Padua points in integrating the Dixon–Szegö function:

In[7]:=
ResourceFunction[
 "PaduaIntegrate"][(4 - 2.1 x^2 + x^4/3) x^2 + x y + 4 (y^2 - 1) y^2, {x, -2, 2}, {y, -1.25, 1.25}, "PaduaType" -> 3]
Out[7]=

WorkingPrecision (1) 

Use 25-digit precision for the numerical integration:

In[8]:=
ResourceFunction[
 "PaduaIntegrate"][(4 - 21/10 x^2 + x^4/3) x^2 + x y + 4 (y^2 - 1) y^2, {x, -2, 2}, {y, -(5/4), 5/4}, WorkingPrecision -> 25]
Out[8]=

Properties and Relations (2) 

If the input function is a polynomial of degree k, PaduaIntegrate gives the exact answer provided k is less than or equal to the setting for InterpolationOrder:

In[9]:=
ResourceFunction[
 "PaduaIntegrate"][-(x^2 y - x - 1)^2 - (x^2 - 1)^2, {x, -3, 3}, {y, -3, 3}, InterpolationOrder -> 6]
Out[9]=
In[10]:=
% - Integrate[-(x^2 y - x - 1)^2 - (x^2 - 1)^2, {x, -3, 3}, {y, -3, 3}]
Out[10]=

A test function due to Franke:

In[11]:=
franke[x_, y_] := 3/4 Exp[-(((9 x - 2)^2 + (9 y - 2)^2)/4)] + 3/4 Exp[-((9 x + 1)^2/49) - (9 y + 1)/10] + 1/2 Exp[-(((9 x - 7)^2 + (9 y - 3)^2)/4)] - 1/5 Exp[-(9 x - 4)^2 - (9 y - 7)^2]

Integrate the function by integrating the interpolant obtained from the resource function PaduaInterpolation:

In[12]:=
pin = ResourceFunction["PaduaInterpolation"][
   franke[x, y], {x, 0, 1}, {y, 0, 1}, InterpolationOrder -> 19, WorkingPrecision -> 30];
Integrate[Integrate[pin[x, y], {x, 0, 1}], {y, 0, 1}]
Out[13]=

The result of directly using PaduaIntegrate is more accurate:

In[14]:=
ResourceFunction["PaduaIntegrate"][franke[x, y], {x, 0, 1}, {y, 0, 1},
  InterpolationOrder -> 19, WorkingPrecision -> 30]
Out[14]=

Version History

  • 1.0.0 – 22 April 2021

Source Metadata

Related Resources

License Information