# Numerical evidence of Eq. (2.25)

The Jacobi polynomials satisfy the recursion relation

\begin{equation}
    %(1-x^2) \partial^k \boldsymbol{P} &= H_{[k]} B H_{[k-1]}^{+} \partial^{k-1} \boldsymbol{P}, \label{eq:recurB0}\\
    (1-x^2)^k \partial^k \boldsymbol{P} = (C^{(k)})^T  \boldsymbol{P}, \label{eq:recurB1}
\end{equation}

where $(C^{(k)})^T = H^{(k)} B^k H^{-1}$ and $H^{(k)} = \text{diags}(h^{(k,\alpha,\beta)}_0, h^{(k,\alpha,\beta)}_1, \ldots)$, $H=H^{(0)}$, and $h^{(k,\alpha,\beta)}_n$ is given as

\begin{equation}
    \left( \partial^k P^{(\alpha, \beta)}_{n},  \partial^k P^{(\alpha, \beta)}_{n}\right)_{\omega^{(\alpha+k,\beta+k)}} = h^{(k,\alpha, \beta)}_{n}, \quad \text{for} \quad n\ge 0.
    \label{eq:orthoderP}
\end{equation}

On index form we have

\begin{equation}
    (1-x^2)^k \partial^k {P}^{(\alpha,\beta)}_n = \sum_{m=n-k}^{n+k}\frac{h^{(k,\alpha, \beta)}_{n}}{h^{(0,\alpha, \beta)}_{m}} b^{(k)}_{nm}{P}_{m}^{(\alpha,\beta)}, \label{eq:recurB2}
\end{equation}

We can shown that this relations holds for some random choices of parameters. If the code cell below runs without error, then the check passes for Legendre and Chebyshev polynomials.

# Implementation


```python
import numpy as np
import sympy as sp
from shenfun.jacobi.recursions import alfa, beta, qn, dqn, n, cn, un, x, h, b, matpow, half

def lhs(n, k, alfa, beta, gn=1):
    """Left hand side of (1)
    
    Parameters
    ----------
    n : int
        Index
    k : int
        parameter
    alfa, beta : numbers
        Jacobi parameters
    gn : scaling function
    """
    return (1-x**2)**k*dqn(alfa, beta, n, k, gn)

def rhs(n, k, alfa, beta, gn=1):
    """Right hand side of (1)
    
    Matrix vector product for given row n.
    
    Parameters
    ----------
    n : int
        Index
    k : int
        parameter
    alfa, beta : numbers
        Jacobi parameters
    gn : scaling function
    """
    s = 0
    for j in range(-k, k+1):
        m = n+j
        s += h(alfa, beta, n, k, gn)/h(alfa, beta, m, 0, gn)*matpow(b, k, alfa, beta, n, m, gn)*qn(alfa, beta, m, gn)
    return s

# Sympy is not smart enough to simplify lhs and rhs to the same expression. Hence
# we need to show numerically that the expressions are zero for some random x and n.

def assert_equal(n, k, alf, bet, x_=1, gn=1):
    f = abs(lhs(n, k, alf, bet, gn) - rhs(n, k, alf, bet, gn)).subs(x, x_)
    assert f < 1e-9
    
# Legendre polynomials with gn=1
for k in (1, 2, 3, 4):
    for x_ in np.random.randn(2):
        assert_equal(6, k, 0, 0, x_)

# Chebyshev polynomials with gn=cn
for k in (1, 2, 3, 4):
    for x_ in np.random.randn(2):
        assert_equal(6, k, -half, -half, x_, cn)

        
```

8aa94ae5794460a62f598bfe2310fc13981c9da6
ipynb
Jupyter Notebook
binder/Equation (2.25).ipynb
spectralDNS/PG-paper-2022
0bfa82e1ff77a5bb9a6ec5be930b4c00b449a1e0
BSD-2-Clause
binder/Equation (2.25).ipynb
spectralDNS/PG-paper-2022
0bfa82e1ff77a5bb9a6ec5be930b4c00b449a1e0
BSD-2-Clause
binder/Equation (2.25).ipynb
spectralDNS/PG-paper-2022
0bfa82e1ff77a5bb9a6ec5be930b4c00b449a1e0
BSD-2-Clause
Qwen/Qwen-72B
1. YES
2. YES