Jacobi Orthogonal Polynomials - Scipy (Python) vs My Recurrence Implementation

Question

I implemented the Jacobi orthogonal polynomial recurrence relation using the following:

Source: http://lsec.cc.ac.cn/~hyu/teaching/shonm2013/STWchap3.2p.pdf

However, when I compare to the Scipy (Python package) implementation, the 2 results are very slightly different for the same inputs:

I have double-checked my recurrence relation and it matches the formula source. How can I determine which result is "correct"? And why they would be slightly different?

Code provided below for anyone Python inclined:

import numpy as np
from scipy.special import jacobi
import matplotlib.pyplot as plt
import matplotlib
import seaborn as sns


cmap = matplotlib.colormaps["tab10"]
#cmap = plt.cm.get_cmap('tab10')
colors = [cmap(i) for i in range(cmap.N)]

#%% Jacobi Polynomials
x = np.arange(-1.0, 1.0, 0.01)
#x = 0.25
alpha = 1
beta = 1
power = 5 # max polynomial power to derive (order)

degrees = np.arange(0,power+1,1)

sns.set(rc = {"figure.dpi":300, 'savefig.dpi':300, 'figure.figsize':(8,6)})
sns.set_theme()
sns.set_style("whitegrid",{'axes.linewidth': 2, 'axes.edgecolor':'black'})

fig,ax = plt.subplots(nrows=1, ncols=1)

# ---- Scipy Implementation
J_scipy = []
for i,degree in enumerate(degrees):
    P = jacobi(degree,alpha,beta,monic=False)(x)
    J_scipy.append(P)
    ax.plot(x,P,label=f"J_scipy_{i}",color=colors[i])

# ---- My recurrence implementation
J0 = np.ones(len(x))
J1 = 0.5*(alpha + beta + 2)*x + 0.5*(alpha - beta)
ax.plot(x,J0,label="J_0",color=colors[0],linestyle="--")
ax.plot(x,J1,label="J_1",color=colors[1],linestyle="--")
J = [J0, J1] # start list of J's
for n in range(2,power+1):
    an = (2*n+alpha+beta+1)*(2*n+alpha+beta+2)/(2*(n+1)*(n+alpha+beta+1))
    bn = (beta**2 - alpha**2)*(2*n+alpha+beta+1)/(2*(n+1)*(n+alpha+beta+1)*(2*n+alpha+beta))
    cn = (n+alpha)*(n+beta)*(2*n+alpha+beta+2)/((n+1)*(n+alpha+beta+1)*(2*n+alpha+beta))
    
    Jn = (an*x - bn)*J[n-1] - cn*J[n-2]
    J.append(Jn)
    
    ax.plot(x,Jn,label=f"J_{n}",color=colors[n],linestyle="--")
    
plt.legend()
plt.title("Jacobi Polynomials")
plt.xlabel("X")
plt.ylabel("P_n")

plt.show(); plt.close()

Answer 1

The recurrence relation is implemented wrongly. The one you show in the document has the following form

$$ J_{n+1} = \left(a_n^{\alpha, \beta}x - b_n^{\alpha,\beta}\right) J_n - c_n^{\alpha,\beta} J_{n-1} $$

In you implementation, you are writing $$ J_{n} = \left(a_n^{\alpha, \beta}x - b_n^{\alpha,\beta}\right) J_{n-1} - c_n^{\alpha,\beta} J_{n-2} $$

The degrees for the polynomial changed, but the $n$ in the coefficients were not changed accordingly.

You might get some differences between different implementations of the recurrence relation, but the error from that would be close to machine precision, typically $O(10^{-13})$ to $O(10^{-14})$. Such errors should not be visible when you are inspecting the curves at this scale.

Answer 2

It turns out there are slightly different derivations/implementations of the recurrence relation. I have identified 3 different implementations so far. Scipy follows the implementation shown here: https://mathworld.wolfram.com/JacobiPolynomial.html

So the difference was due to different recurrence relation derivations.