Repository: scipy/scipy
File: scipy/special/_orthogonal.py (Python)

Code before:

    >>> import numpy as np
    >>> from scipy.special import sh_legendre, legendre
    >>> x = np.arange(0.0, 1.0, 0.01)
    >>> n = 3
    >>> np.allclose(sh_legendre(n)(x), legendre(n)(2*x - 1))
    True

    The polynomials :math:`P_n^*` satisfy a recurrence relation obtained by the change of variables :math:`t = 2x - 1` in the standard Legendre recurrence:
    
    .. math::
    
        (n+1) P_{n+1}^*(x) = (2n+1)(2x-1)\,P_n^*(x) - n\,P_{n-1}^*(x).

    This can be easily checked on :math:`[0, 1]` for :math:`n = 3`:
    >>> x = np.arange(0.0, 1.0, 0.01)
    >>> lhs = (3 + 1) * sh_legendre(4)(x)
    >>> rhs = (2*3 + 1) * (2*x - 1) * sh_legendre(3)(x) - 3 * sh_legendre(2)(x)
    >>> np.allclose(lhs, rhs)
    True

    Orthogonality over :math:`[0,1]` with weight 1 can be checked numerically; for example, :math:`P_2^*` is orthogonal to :math:`P_3^*`:

    >>> x = np.linspace(0.0, 1.0, 400)
    >>> y = sh_legendre(2)(x) * sh_legendre(3)(x)
    >>> np.isclose(np.trapz(y, x), 0.0, atol=1e-4)
    True

    See Also
    --------
    scipy.special._orthogonal.legendre
    scipy.special._orthogonal.roots_sh_legendre
    """
    if n < 0:
        raise ValueError("n must be nonnegative.")

    def wfunc(x):
        return 0.0 * x + 1.0
    if n == 0:
        return orthopoly1d([], [], 1.0, 1.0, wfunc, (0, 1), monic,
                           lambda x: _ufuncs.eval_sh_legendre(n, x))
    x, w = roots_sh_legendre(n)
    hn = 1.0 / (2 * n + 1.0)
    kn = _gam(2 * n + 1) / _gam(n + 1)**2
    p = orthopoly1d(x, w, hn, kn, wfunc, limits=(0, 1), monic=monic,
                    eval_func=lambda x: _ufuncs.eval_sh_legendre(n, x))
    return p


# Make the old root function names an alias for the new ones
_modattrs = globals()

Review comment (bug):
Might want to consider using `scipy.integrate.trapezoid()` here. The reasons are not super interesting, but if you're morbidly curious there were a few older complaints:

- https://github.com/numpy/numpy/issues/25586
- https://github.com/numpy/numpy/pull/25738
- https://github.com/scipy/scipy/issues/20486
- https://github.com/scipy/scipy/issues/12924

Some of the CI failures I see here are unrelated to your changes and involve the latest release of Cython `3.2.0` on Nov 5th as described at gh-23932.

I don't know enough about Legendre polynomials to have an opinion on the content beyond that.

Code after:
    to the non-shifted polynomials :math:`P_n` by
    :math:`P_n^*(x) = P_n(2x - 1)`. We can verify this on the
    interval :math:`[0, 1]`:

    >>> import numpy as np
    >>> from scipy.special import sh_legendre, legendre
    >>> from scipy.integrate import trapezoid
    >>> x = np.arange(0.0, 1.0, 0.01)
    >>> n = 3
    >>> np.allclose(sh_legendre(n)(x), legendre(n)(2*x - 1))
    True

    The polynomials :math:`P_n^*` satisfy a recurrence
    relation obtained by the change of variables
    :math:`t = 2x - 1` in the standard Legendre recurrence:
    
    .. math::
    
        (n+1) P_{n+1}^*(x) = (2n+1)(2x-1)\,P_n^*(x) - n\,P_{n-1}^*(x).

    This can be easily checked on :math:`[0, 1]`
    for :math:`n = 3`:

    >>> n = 3
    >>> x = np.linspace(0.0, 1.0, 101)
    >>> lhs = (n + 1) * sh_legendre(n + 1)(x)
    >>> rhs = (
    ...     (2*n + 1) * (2*x - 1) * sh_legendre(n)(x)
    ...     - n * sh_legendre(n - 1)(x)
    ... )
    >>> np.allclose(lhs, rhs)
    True

    Orthogonality over :math:`[0,1]` with weight 1 can be
    checked numerically; for example, :math:`P_2^*`
    is orthogonal to :math:`P_3^*`:

    >>> x = np.linspace(0.0, 1.0, 400)
    >>> y = sh_legendre(2)(x) * sh_legendre(3)(x)
    >>> np.isclose(trapezoid(y, x), 0.0, atol=1e-12)
    True

    See Also
    --------
    scipy.special.legendre
    scipy.special.roots_sh_legendre
    """
    if n < 0:
        raise ValueError("n must be nonnegative.")

    def wfunc(x):