A library that implements Jacobi polynomials and Gauss quadrature related operations.
In the notebooks directory there are tow notebooks, as follows, that shows the basic usage of the library.
 Calculation of some orthogonal polynomials and their zeros
 Gausstype qudrature, interpolation and differentiation
Jacobi polynomials jacobi(x, n, a, b)
where x
is where you want to calculate the
polynomial of degree n
with weights a
and b
.
To calculate the derivative of jacobi polynomials, use function
djacobi(x, n, a, b)
.
When calculating the weights of Gaussian qudrature, it is necessary to determine
the zeros of Jacobi polynomials. The function jacobi_zeros(m, a, b)
calculates the
zeros of a Jacobi polynomial of degree m
with weights a
and b
. Also present
is the non allocating function jacobi_zeros!(m, a, b, x)
.
Legendre polynomials are a special case when a and b are zero.
Function legendre(x, n)
implements this simpler recurrence relation.
Also available are Chebyshev polynomials and the respective derivatives and zeros:

chebyshev(x, n)
Chebyshev polynomial of the first kind$T_n(x)$ . 
chebyshev2(x, n)
Chebyshev polynomial of the second kind$U_n(x)$ . 
chebyshev_zeros(n)
andchebyshev2_zeros(n)
Roots of Chebyshev polynomials 
chebyshev_zeros!(n, x)
andchebyshev2_zeros!(n, x)
Roots of Chebyshev polynomials
Gaussian quadrature is commonly used to evaluate integrals numerically. To evaluate integrals using Gaussian quadrature, a set of nodes and its corresponding weights must be known. In the most basic algorithm, no node can be specified by the user. It is often convenient to include one or both ends of the domain and this provides different quadrature rules.
The different quadrature rules implemented are:
 GaussJacobi, where no nodes are specified (extension
gj
).  GaussLobattoJacobi, where both ends of the domain are included (extension
glj
).  GaussRadauJacobi, where the node 1 (leftmost end of the domain) is included
(extension
grjm
).  GaussRadauJacobi, where the node +1 (rightmost end of the domain) is included
(extension
grjp
).
These formulations implement the numerical integration of $$ \int_{1}^1 (1x)^a(1+x)^b f(x)dx $$
The functions zgj
, zglj
, zgrjm
and zgrjp
calculate the nodes of the quadrature
rule. The corresponding weights are calculated using functions wgj
, wglj
,
wgrjm
and wgrjp
. Of course in both cases, the weights a
and b
must be specified.
As an example
fun(x) = sin(x)
z = zglj(5, 0.0, 0.0)
w = wglj(z, 0.0, 0.0)
f = fun(z)
Ix = sum(w .* f)
where Ix
is the estimate of the integral.
Gaussian quadrature is useful because it allows the exact integration of polynomials using few nodes (the exact order depends on the quadrature rule cited above).
But there is another advantage: the nodes of the quadrature rule is very convenient for interpolating functions using high order polynomials and is commonly used in high order finite element procedures such as hpFEM or spectral element method. In these applications, it is necessary to compute derivatives and interpolate data from different grids.
The functions dgj
, dglj
, dgrjm
and drjp
calculates the derivative matrix as
shown in the example that follows:
fun(x) = sin(x)
z = zglj(5, 0.0, 0.0)
D = dglj(z, 0.0, 0.0)
f = fun(z)
df = D * f
where df
is an estimate of the derivative of the function at the quadrature nodes.
Another important operation is interpolation. If a function is known at some nodes,
in this case the quadrature nodes, how can we accurately interpolate the function on
other nodes? Since we know the nodes, Lagrangian interpolation is the best way.
The function lagrange
implements the standard definition of the Langrangian
interpolation. The example below plots the Lagrangian interpolators of
the GaussLobattoJacobi quadrature points for 5 nodes.
using PyPlot
Q = 5
z = zglj(Q)
nx = 201
x = 1.0:0.01:1.0
y = zeros(nx, Q)
for k = 1:Q, i=1:nx
y[i,k] = lagrange(k, x[i], z)
end
for k=1:Q
plot(x, y[:,k])
end
If the operation above is to be repeated often, precalculating the Lagrangian
interpolators is useful and an Interpolation matrix can be calculated. The
following example illustrates the use of the interpolation matrix that can
be computed with the function interp_mat
.
using PyPlot
Q = 5
z = zglj(Q)
nx = 201
x = 1.0:0.01:1.0
ye = sin(pi*z)
ye2 = sin(pi*x)
Im = interp_mat(x, z)
y = Im * ye
plot(z, ye, "o")
plot(x, ye2, "r")
plot(x, y, "b")
increasing the number of quadrature points the interpolated function (blue line) becomes more accurate.
This package was implemented using both references below.
 Spectral/hp Element Methods for CFD, 2nd edition, Karniadakis and Sherwin, 2005.
 NIST Handbook of Mathematical Functions (http://dlmf.nist.gov/18)