Legendre Polynomial

legendre.py

"""
Computes Legendre Polynomials


         M      m      (2n-2m)!        n-2m
P (x) = sum (-1)  ------------------- x
 n      m=0        n
                  2 m! (n-m)! (n-2m)!

where M = n // 2

"""

import sympy
from sympy import Rational, factorial, symbols, exp
from pprint import pprint
import matplotlib.pyplot as plt


def legendre_polynomial(n):
    M = n // 2

    x, P = symbols('x P')
    for m in range(0, M+1):
        term = (-1)**m * factorial(2*n-2*m) / (2**n*factorial(m)*factorial(n-m)*factorial(n-2*m))*x**(n-2*m)
                                       
        if P == symbols('P'):
            P = term
        else:
            P += term

    return P, x


def plot_curve(y, x, up=1, low=-1, num_points=50, label="curve"):
    """
    known_coeffs (list): list of tuples (symbol, value) that is used
        to substitute coefficients in `y`.
    """
    # Plot, between -2 to 2, the value of the series for up to x^7.
    step_size = (up - low) / num_points
    xvals = []
    yvals = []

    xval = low
    while xval < up:
        yval = float(y.subs([(x, xval)]))
        xvals.append(xval)
        yvals.append(yval)
        xval += step_size

    plt.xlabel('x')
    plt.ylabel('y')

    plt.plot(xvals, yvals, "o-", label=label)



if __name__ == "__main__":
    P6, x1 = legendre_polynomial(6)
    P8, x2 = legendre_polynomial(8)

    print("P6 = " + str(P6))
    print("P8 = " + str(P8))

    plot_curve(P6, x1, label="$P_6$")
    plot_curve(P8, x2, label="$P_6$")
    plt.legend(loc="upper right")
    plt.show()