Source code for Gauss_Lobatto_nodes

# -*- coding: utf-8 -*-
"""⭕ To access the source code, click on the [source] button at the right
side or click on
:download:`[Gauss_Lobatto_nodes.py]</contents/TEACHING/advanced_numerical_methods/mimetic_spectral_element_method/Gauss_Lobatto_nodes.py>`.
Dependence may exist. In case of error, check import and install required
packages or download required scripts. © mathischeap.com
"""

import numpy as np
from scipy.special import legendre


def _newton_method_(f, df_dx, x_0, n_max):
    """Newton method for root finding."""
    x = [x_0,]
    for i in range(n_max - 1):
        # noinspection PyUnresolvedReferences
        x.append(x[i] - f(x[i]) / df_dx(x[i]))
        if abs(x[i + 1] - x[i]) < np.finfo(float).eps * 10:
            break
    return x[-1]




def Gauss_Lobatto_nodes(p):
    """Compute the Gauss Lobatto nodes of degree ``p``
    using newton method. Credits belong to Lorenzo.

    :param p: The quadrature degree (positive integer).
    :type p: int
    :return:

        #. (numpy.ndarray) - The Lobatto quadrature nodes.

    :example:

    >>> n = Gauss_Lobatto_nodes(3)
    >>> n
    array([-1.       , -0.4472136,  0.4472136,  1.       ])
    """
    x_0 = np.cos(np.arange(1, p) / p * np.pi)
    nodal_pts = np.zeros((p + 1))
    nodal_pts[0] = 1
    nodal_pts[-1] = -1
    legendre_prime = lambda x, n: (n * legendre(n - 1)(x) - n * x * legendre(n)(x)) / (1 - x ** 2)
    for i, ch_pt in enumerate(x_0):
        leg_p = lambda x: (1 - x**2)**2 * legendre_prime(x,p)
        leg_pp = lambda x: (2 * x * legendre_prime(x, p) - p * (p + 1) * legendre(p)(x)) * (1 - x ** 2)
        nodal_pts[i + 1] = _newton_method_(leg_p, leg_pp, ch_pt, 100)

    return nodal_pts[::-1]



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