"""
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
packages or download required scripts. © mathischeap.com
"""
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()
