In [1]:
import numpy as np
import matplotlib.pyplot as plt
import numpy.linalg
In [2]:
np.set_printoptions(precision=3)
In [3]:
# How many Legendre polynomials do you want?
n = 6

Matrix representations¶

Matrix representation of a polynomial and its derivatives

$$ P(x) = a_0 + a_1 x + a_2 x^2 + ... + a_n x^n \to (a_0, a_1, a_2, ..., a_n)^T$$

$$ \frac{d}{dx}P(x) = a_1 + 2a_2x + 3a_3x^2 + ... + n a_n x^{n-1} \to \left( \begin{matrix} 0 & 1 & 0 &...& 0 \\ 0 & 0 & 2 &...& 0 \\ \vdots & & & & \\ 0 & 0 & 0 &...& n \\ 0 & 0 & 0 &...& 0 \end{matrix} \right)$$

$$ \frac{d^2}{dx^2}P(x) = 2a_2 + 6a_3x + ... + n (n-1) a_n x^{n-2} \to \left( \begin{matrix} 0 & 0 & 2 & 0 &...& 0 \\ 0 & 0 & 0 & 6 &...& 0 \\ \vdots & & & & & \\ 0 & 0 & 0 & 0 &...& n(n-1) \\ 0 & 0 & 0 & 0 &...& 0 \\ 0 & 0 & 0 & 0 &...& 0 \end{matrix} \right)$$

Effects of multiplying $x$ and $x^2$ on the derivatives

$$ x\frac{d}{dx}P(x) = a_1x + 2a_2x^2 + 3a_3x^3 + ... + n a_n x^{n} \to \left( \begin{matrix} 0 & 0 & 0 &...& 0 \\ 0 & 1 & 0 &...& 0 \\ 0 & 0 & 2 &...& 0 \\ \vdots & & & & \\ 0 & 0 & 0 &...& n \end{matrix} \right)$$

$$ x^2 \frac{d^2}{dx^2}P(x) = 2a_2 x^2 + 6a_3x^3 + ... + n (n-1) a_n x^{n} \to \left( \begin{matrix} 0 & 0 & 0 & 0 &...& 0 \\ 0 & 0 & 0 & 0 &...& 0 \\ 0 & 0 & 2 & 0 &...& 0 \\ 0 & 0 & 0 & 6 &...& 0 \\ \vdots & & & & & \\ 0 & 0 & 0 & 0 &...& n(n-1) \\ \end{matrix} \right)$$

How to code these matrix represetations? Let's see what the function np.eye does

In [4]:
np.eye(n)
Out[4]:
array([[1., 0., 0., 0., 0., 0.],
       [0., 1., 0., 0., 0., 0.],
       [0., 0., 1., 0., 0., 0.],
       [0., 0., 0., 1., 0., 0.],
       [0., 0., 0., 0., 1., 0.],
       [0., 0., 0., 0., 0., 1.]])
In [5]:
np.eye(n,k=1)
Out[5]:
array([[0., 1., 0., 0., 0., 0.],
       [0., 0., 1., 0., 0., 0.],
       [0., 0., 0., 1., 0., 0.],
       [0., 0., 0., 0., 1., 0.],
       [0., 0., 0., 0., 0., 1.],
       [0., 0., 0., 0., 0., 0.]])

Legendre operator as a matrix¶

$$ \mathrm{Legendre} = (1-x^2) \frac{d^2}{dx^2} - 2x \frac{d}{dx} $$

In [6]:
op1 = np.arange(0,n)*np.eye(n) # x * d/dx
In [7]:
op2 = (np.arange(-1,n-1)*np.arange(0,n))*np.eye(n,k=2)  # d^2/dx^2
In [8]:
op3 = (np.arange(-1,n-1)*np.arange(0,n))*np.eye(n)  # x^2 * d^2/dx^2

Building the full operator

In [9]:
Legendre = op2 - op3 - 2*op1
In [10]:
Legendre
Out[10]:
array([[  0.,   0.,   2.,   0.,   0.,   0.],
       [  0.,  -2.,   0.,   6.,   0.,   0.],
       [  0.,   0.,  -6.,   0.,  12.,   0.],
       [  0.,   0.,   0., -12.,   0.,  20.],
       [  0.,   0.,   0.,   0., -20.,   0.],
       [  0.,   0.,   0.,   0.,   0., -30.]])

Computation of eigenvalues and eigenvectors¶

In [11]:
eigval, eigvec = numpy.linalg.eig(Legendre)
In [12]:
eigval
Out[12]:
array([  0.,  -2.,  -6., -12., -20., -30.])

Each column of the eigenvectors gives the coefficients of the Legendre polynomials (not normalized)

In [13]:
eigvec
Out[13]:
array([[ 1.   ,  0.   , -0.316,  0.   ,  0.065,  0.   ],
       [ 0.   ,  1.   ,  0.   , -0.514,  0.   ,  0.157],
       [ 0.   ,  0.   ,  0.949,  0.   , -0.649,  0.   ],
       [ 0.   ,  0.   ,  0.   ,  0.857,  0.   , -0.734],
       [ 0.   ,  0.   ,  0.   ,  0.   ,  0.758,  0.   ],
       [ 0.   ,  0.   ,  0.   ,  0.   ,  0.   ,  0.661]])

Discretization $x$ for plotting the polynomials

In [14]:
nx = 101
x = np.linspace(-1,1,nx)
In [15]:
exponent = np.arange(0,n)

Building the set of monomials $\{0, x, x^2, x^3, ...\}$

In [16]:
monomials = (x[:,np.newaxis]**exponent)

Building the polynomial eigenfunctions

In [17]:
polynomials = monomials @ (eigvec)

Plotting the eigenfunctions (Legendre polynomials)

In [18]:
for i in range(n):
    plt.plot(x,polynomials[:,i]/np.abs(polynomials[0,i]),label=rf"$P_{i}(x)$")
plt.legend()
Out[18]:
<matplotlib.legend.Legend at 0x109a6fcb0>
No description has been provided for this image