Function Repository Resource:

PaduaPoints

Source Notebook

Generate Padua points for bivariate interpolation and cubature

Contributed by: Jan Mangaldan

ResourceFunction["PaduaPoints"][n]

gives the type 1 Padua points of degree n over the domain [-1,1]×[-1,1].

ResourceFunction["PaduaPoints"][m,n]

gives the type m Padua points over the domain [-1,1]×[-1,1].

ResourceFunction["PaduaPoints"][m,n,{{xmin,xmax},{ymin,ymax}}]

gives the Padua points over the domain [xmin,xmax]×[ymin,ymax].

ResourceFunction["PaduaPoints"][m,n,{{xmin,xmax},{ymin,ymax}},prec]

uses the working precision prec.

Details

The Padua points are a unisolvent point set that can be used to interpolate and integrate a bivariate function over a rectangular domain.
Possible types m of Padua points are 1, 2, 3 and 4.
The result of ResourceFunction["PaduaPoints"] is a list of length (n+1)(n+2)/2.

Examples

Basic Examples (3) 

Padua points of degree 3:

In[1]:=
ResourceFunction["PaduaPoints"][3]
Out[1]=

Padua points of type 4 and degree 3:

In[2]:=
ResourceFunction["PaduaPoints"][4, 3]
Out[2]=

Padua points of type 1 and degree 3 over[0,1]×[0,1]:

In[3]:=
ResourceFunction["PaduaPoints"][1, 3, {{0, 1}, {0, 1}}]
Out[3]=

Use 25-digit precision:

In[4]:=
ResourceFunction["PaduaPoints"][1, 3, {{0, 1}, {0, 1}}, 25]
Out[4]=

Applications (2) 

A function to approximate:

In[5]:=
peaks[x_, y_] := 3 (1 - x)^2 E^(-x^2 - (y + 1)^2) - 10 (x/5 - x^3 - y^5) E^(-x^2 - y^2) - 1/3 E^(-(x + 1)^2 - y^2)

Construct a degree 20 Padua interpolant over[-3,3]×[-3,3]:

In[6]:=
padua = ResourceFunction["PaduaPoints"][1, 20, {{-3, 3}, {-3, 3}}, 25];
approx[x_, y_] = InterpolatingPolynomial[
   Transpose[{padua, peaks @@@ padua}], {x, y}];

Compare the interpolant with the original function:

In[7]:=
{Plot3D[approx[x, y], {x, -3, 3}, {y, -3, 3}, PlotRange -> All], Plot3D[peaks[x, y], {x, -3, 3}, {y, -3, 3}, PlotRange -> All]} // GraphicsRow
Out[7]=

Plot the difference:

In[8]:=
Plot3D[approx[x, y] - peaks[x, y], {x, -3, 3}, {y, -3, 3}, PlotRange -> All]
Out[8]=

Approximate the integral of a Gaussian function using a degree-25 Padua approximation:

In[9]:=
gau = {x, y} |-> Exp[-x^2 - y^2];
padua = ResourceFunction["PaduaPoints"][1, 25, {{0, 1}, {0, 1}}];
Integrate[
 Integrate[
  InterpolatingPolynomial[
   Transpose[{padua, gau @@@ padua}], {x, y}], {y, 0, 1}], {x, 0, 1}]
Out[9]=

Compare with the exact result:

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

Properties and Relations (1) 

The Padua points lie on a Lissajous curve:

In[11]:=
With[{n = 7},
 ParametricPlot[{-Cos[(n + 1) t], -Cos[n t]}, {t, 0, \[Pi]}, Epilog -> {Directive[ColorData[97, 4], AbsolutePointSize[4]], Point[ResourceFunction["PaduaPoints"][n]]}]]
Out[11]=

Version History

  • 1.0.0 – 05 April 2021

Source Metadata

Related Resources

License Information