# Source code for Lagrange_and_edge_polynomials

"""
Lagrange polynomials and edge polynomials in  the 1d reference element
:math:\Omega_{\mathrm{ref}}=[-1,1].

⭕ To access the source code, click on the [source] button at the right
side or click on
:download:[Lagrange_and_edge_polynomials.py]</contents/LIBRARY/ptc/mathischeap_ptc/Lagrange_and_edge_polynomials.py>.
Dependence may exist. In case of error, check import and install required
"""

import numpy as np

[docs]
def Lagrange_polynomials(nodes, xi):
"""Compute the Lagrange polynomials :math:l^{i}(\\xi) built on
nodes. They are then evaluated at xi.

:param nodes: The nodes on which the Lagrange polynomials are built.
:type nodes: 1d np.array
:param xi: The points to evaluate the Lagrange polynomials.
:type xi: 1d np.array
:return: A 2d numpy.array whose 0-dimension refers to the indexes of
the polynomials and 1-dimension refers to the values evaluated
at xi.
:example:

>>> nodes = np.array([-1., -0.4472136,  0.4472136,  1.])
>>> x = np.linspace(-1, 1, 50)
>>> lp = Lagrange_polynomials(nodes, x)
>>> lp[0,:] # l^0(xi) # doctest: +ELLIPSIS
array([ 1.        ,  0.88167345,  0.77142177,  0.66898996,...
"""
p = np.size(nodes)
basis = np.ones((p, np.size(xi)))
for i in range(p):
for j in range(p):
if i != j:
basis[i, :] *= (xi - nodes[j]) / (nodes[i] - nodes[j])
return basis

[docs]
def edge_polynomials(nodes, xi):
"""Compute the edge polynomials :math:e^{i}(\\xi) built on
nodes. They are then evaluated at xi.

:param nodes: The nodes on which the edge polynomials are built.
:type nodes: 1d np.array
:param xi: The points to evaluate the edge polynomials.
:type xi: 1d np.array
:return: A 2d numpy.array whose 0-dimension refers to the indexes of
the polynomials and 1-dimension refers to the values evaluated
at xi.
:example:

>>> nodes = np.array([-1., -0.65465367,  0.,  0.65465367,  1.])
>>> x = np.linspace(-1, 1, 50)
>>> ep = edge_polynomials(nodes, x)
>>> ep[2,:]  # e^3(xi) # doctest: +ELLIPSIS
array([ 0.91016418,  0.60853048,  0.34702611,  0.12374712, ...
"""
p = np.size(nodes) - 1
xi_xj = nodes.reshape(p + 1, 1) - nodes.reshape(1, p + 1)
xi_xj[np.diag_indices(p + 1)] = 1
c_i = np.prod(xi_xj, axis=1)
c_i_div_cj = np.transpose(c_i.reshape(1, p+1) / c_i.reshape(p+1, 1))
derivative = c_i_div_cj / xi_xj
derivative[np.diag_indices(p + 1)] = 0
derivative[np.diag_indices(p + 1)] = -np.sum(derivative, axis=1)
polynomials = Lagrange_polynomials(nodes, xi)
derivatives_poly = np.transpose(derivative) @ polynomials
edge_poly = np.zeros((p, np.size(xi)))
for i in range(p):
for j in range(i + 1):
edge_poly[i] -= derivatives_poly[j, :]
return edge_poly

if __name__ == '__main__':
import doctest
doctest.testmod()

import matplotlib.pyplot as plt

nodes = np.array([-1., -0.65465367, 0., 0.65465367, 1.])
x = np.linspace(-1, 1, 200)
lp = Lagrange_polynomials(nodes, x)
ep = edge_polynomials(nodes, x)

fig = plt.figure(figsize=(16,6))

ax = plt.subplot2grid((1, 2), (0, 0))
for i in range(5):
plt.plot(x, lp[i,:])
plt.title('Lagrange polynomials')

ax = plt.subplot2grid((1, 2), (0, 1))
for i in range(4):
plt.plot(x, ep[i,:])
plt.title('Edge polynomials')

plt.show()