Preface


This section presents the basic properties of solutions to the Laplace and Helmholtz equations.

Return to computing page for the first course APMA0330
Return to computing page for the second course APMA0340
Return to Mathematica tutorial for the first course APMA0330
Return to Mathematica tutorial for the second course APMA0340
Return to the main page for the first course APMA0330
Return to the main page for the second course APMA0340
Return to Part V of the course APMA0340
Introduction to Linear Algebra with Mathematica

Legendre Polynomials


Adrien-Marie Legendre

Adrien-Marie Legendre (1752--1833) was a French mathematician. Legendre made numerous contributions to mathematics. Well-known and important concepts such as the Legendre polynomials and Legendre transformation are named after him.

Consider Legendre's equation

\[ \left( 1-x^2\right) y'' -2x\,y' + n(n+1)\, y =0 , \]
which can be rewritten in self-adjoint form
\[ \frac{\text d}{{\text d}x} \left[ \left( 1- x^2 \right) \frac{{\text d} y}{{\text d}x} \right] + n(n+1)\, y =0 , \]
where n is a real constant, but it is usually an integer. The above equation is named after a French mathematician Adrien-Marie Legendre (1752--1833) who introduced the Legendre polynomials in 1782. Legendre's equation comes up in many physical situations involving spherical symmetry. When n is an integer, the Legandre differential equation has a polynomial solution (with the normalization \( P_n (1) =1 \) ) that is usually denoted by \( P_n (x) \) and is called Legendre's polynomial of order n. The Legendre equation has non-polynomial solutions represented by infinite series. These are the Legendre functions of the second kind, denoted by \( Q_n (x) . \) Each Legendre polynomial may be expressed using Rodrigues' formula (named after Benjamin Olinde Rodrigues (1795--1851), a French banker, mathematician, and social reformer):
\[ P_n (x) = \frac{1}{2^n n!} \,\frac{{\text d}^n}{{\text d}x^n} \left( x^2 -1 \right)^n , \qquad n=1,2,\ldots . \]
They can also be defined via Bonnet’s recursion formula:
\[ \left( n+1 \right) P_{n+1} (x) = \left( 2n+1 \right) x\, P_n (x) - n\,P_{n-1} (x), \qquad n=2,3, \ldots ; \quad P_0 (x) =1, \quad P_1 (x) =x; \]
which follow from the generating function for Legendre polynomials:
\[ \frac{1}{\sqrt{1- 2xt + t^2}} = \sum_{n\ge 0} P_n (x)\, t^n \qquad\Longrightarrow \qquad P_n (x) = \frac{1}{2^n} \,\sum_{k=0}^n \binom{n}{k}^2 \left( x-1 \right)^{n-k} \left( x+1 \right)^k . \]
Its exponential generating function is
\[ e^{xz}\,J_0 \left( z \sqrt{1-x^2} \right) = \sum_{n\ge 0} P_n (x)\, \frac{z^n}{n!} . \]
The Rodrigues formula can be derived from the Lagrange inversion theorem (LIT) using
\[ w = f(z) = 2 \, \frac{z- z_0}{z^2 -1} , \]
where \( z_0 = z(0) , \) with \( w_0 = f(z_0 ) . \) Then
\[ \frac{1}{1 - (w- w_0 ) \, \frac{{\text d}}{{\text d} x} \left( \frac{z-z_0}{f(z) - w_0} \right)_{z=z(w)}} = \sum_{n\ge 0} c_n \left( w- w_0 \right)^n , \]
where
\[ c_n = \frac{1}{n!} \, \frac{{\text d}^n}{{\text d} z^n} \left. \left( \frac{z- z_0}{f(z) - w_0} \right) \right\vert_{z= z_0} = \frac{1}{2^n \, n!} \, \frac{{\text d}^n}{{\text d} z^n} \left( z^2 -1 \right)^n . \]
Legendre polynomials are eigenfunctions of a singular Sturm–Liouville problem:
\[ \frac{\text d}{{\text d}x} \left[ \left( 1- x^2 \right) \frac{{\text d} y}{{\text d}x} \right] = - \lambda\, y , \qquad y(\pm 1) < \infty , \]
where the eigenvalue λ corresponds to n(n + 1).

Legendre polynomials are symmetric or antisymmetric, that is,

\[ P_n (-x) = (-1)^n P_n (x) , \qquad n=0,1,2,\ldots . \]
We remind the numerical values:
\[ P_n (1) = 1 , \qquad P_n (-1) = (-1)^n , \qquad P_{2n} (0) = (-1)^n \,\frac{(2n)!}{2^{2n} \, (n! )^2} ,\qquad P'_{n} (1) = \frac{n(n+1)}{2} , \qquad n=0,1,2,\ldots . \]
The Legendre polynomials are built into Mathematica. Mathematica's notation is LegendreP[n,x] for \( P_n (x). \) We now use Mathematica to obtain the formulas for the first 11 of these polynomials. We put them in a table.
TableForm[Table[{i, i*(i + 1), LegendreP[i, x]}, {i, 0, 10}], TableHeadings -> {None, {"n", "\[Lambda]_k", "P_n (x)"}}]
To get some idea of what these polynomials look like, we construct graphs of the first 7. We first define a function legraph[n] that produces a graph of the kth polynomial, and then we use a Do loop to construct the first 7 graphs.
legraph[n_] := Plot[LegendreP[n, x], {x, -1, 1}, AxesLabel -> {"x", "P_n (x)"}, PlotRange -> {{-1.1, 1.1}, {-1.1, 1.1}}, PlotLabel -> Row[{"n=", PaddedForm[n, 3]}]]
GraphicsGrid[Table[{legraph[i], legraph[i + 1]}, {i, 0, 6, 2}]]

First we set up a series to work with:

y = Sum[c[i] x^i, {i, 0, 6}] + O[x]^7
Out[1]= c[0] + c[1] x + c[2] x^2 + c[3] x^3 + c[4] x^4 + c[5] x^5 +
c[6] x^6 + O[x]^7

To see that our input is correct we type

OutputForm[y]//Normal

We are going to insert this series into the differential equation

de = (1 -x^2) D[y, {x,2}] - 2 x D[y,x] + 20 y == 0
Out[2]= (20 c[0] + 2 c[2]) + (18 c[1] + 6 c[3]) x + (14 c[2] + 12 c[4])
x^2 + (8 c[3] + 20 c[5]) x^3 + 30 c[6] x^4 + O[x]^5 == 0
Then we get the coefficients
coeffeqns = LogicalExpand[de]
Out[3]= 20 c[0] + 2 c[2] == 0 && 18 c[1] + 6 c[3] == 0 &&
14 c[2] + 12 c[4] == 0 && 8 c[3] + 20 c[5] == 0 && 30 c[6] == 0
and express through first two c[0] and c[1]:
solvedcoeffs = Solve[ coeffeqns, Table[ c[i], {i,2,12}]]

Solve::svars: Equations may not give solutions for all "solve" variables. >>

Out[4]= {{c[2] -> -10 c[0], c[3] -> -3 c[1], c[4] -> (35 c[0])/3, c[5] -> (6 c[1])/5, c[6] -> 0}}
Put these back into the series expansion.
y = y /. solvedcoeffs
Out[5]= { c[0] + c[1] x -10 c[0] x^2 - 3 c[1] x^3 + (35/3) c[0] x^4 + (6/5) c[1] x^5 + O[x]^7 }
Extract out the two linearly independent solutions
Coefficient[ y, c[0]]
Out[6]= {1 - 10 x^2 + (35 x^4)/3}
Coefficient[ y, c[1]]
Out[7]= {x - 3 x^3 + (6 x^5)/5}
To check our answer, we use Mathematica again because it knows Legendre:
LegendreP[4,x]
Out[8]= 1/8 (3 - 30 x^2 + 35 x^4)
This is different from the polynomial obtained earlier, but only by a constant factor. The constant factor is used to set the normailization for the Legendre polynomials.

Legendre Expansions


Since Legendre polynomial are orthogonal
\[ \int_{-1}^1 P_n (x) \, P_m (x) \,{\text d} x = 0 \qquad\mbox{for} \quad n\ne m, \]
any suitable function \( f (x) \) on the interval [-1,1] can be expanded into generalized Fourier--Legendre series:
\[ f(x) = \sum_{n\ge 0} c_n P_n (x) , \]
where the coefficients are
\[ c_n = \left( n + \frac{1}{2} \right) \int_{-1}^1 P_n (x) \, f(x) \,{\text d} x , \qquad n=0,1, \ldots , \]
because
\[ \int_{-1}^1 P_n^2 (x) \,{\text d} x = \frac{2}{2n+1} . \]

Example 1: Find the Fourier-Legendre series expansion of the Heaviside step function H(t), defined on the finite interval [-1,1].

Its series expansion is written in the form

\[ H(t) = \sum_{n\ge 0} c_n \,P_n (t) , \qquad t \in [-1,1] . \]
Substituting the explicit expressions for the Legendre polynomials, we get
\begin{align*} c_n &= \frac{2n+1}{2} \, \int_{-1}^1 H(t)\, P_n (t)\,{\text d}t = \frac{2n+1}{2} \, \int_{0}^1 \, P_n (t)\,{\text d}t . \end{align*}
Integrating by parts, we find the explicit expressions for the coefficients:
\begin{align*} c_n &= \frac{2n+1}{2} \,\frac{1}{2^n \,n!} \int_{0}^1 \frac{{\text d}^n}{{\text d} t^n} \left( t^2 -1 \right)^n {\text d}t = \frac{2n+1}{2^{n+1} n!} \, \left[ \frac{{\text d}^{n-1}}{{\text d} t^{n-1}} \left( t^2 -1 \right)^n \right]_{t=0}^{t=1} . \\ &= \left\{ \begin{array}{ll} 0 , & \ \mbox{if $n$ is even} \\ (-1)^k \,\frac{(4k+3) \,(2k)!}{2^{2k+2}\,k! \,(k+1)!} , & \ \mbox{if $n = 2k+1$ is odd. } \end{array} \right. \end{align*}
Thus, the Fourier--Legendre series expansion of the Heaviside function is given by
\[ H(t) = \frac{1}{2} + \sum_{k\ge 0} \, \frac{(-1)^k \,(4k+3) \,(2k)!}{2^{2k+2}\,k! \,(k+1)!} \, P_{2k+1} (x) . \]
Then we plot finite approximations with 5 and 20 terms:
legendre5 = 1/2 + Sum[(-1)^k *(4*k + 3)* Factorial[2*k]/2^(2*k + 2)/Factorial[k + 1]/Factorial[k] * LegendreP[2*k + 1, t], {k, 0, 5}];
legendre20 = 1/2 + Sum[(-1)^k *(4*k + 3)* Factorial[2*k]/2^(2*k + 2)/Factorial[k + 1]/Factorial[k] * LegendreP[2*k + 1, t], {k, 0, 20}];
Plot[{HeavisideTheta[t], legendre5, legendre20}, {t, -1.1, 1.1}, PlotStyle -> {{Blue, Thick}, Orange, Red}, PlotRange -> {-0.1, 1.2}]

Example 2: Suppose we want to expand the signum function into the Fourier--Legendre series. Using an integration formula

\[ \int_{-1}^x P_n (t) \,{\text d}t = \frac{P_{n+1} (x) - P_{n-1} (x)}{2n+1} , \]
we find the explicit expressions for the coefficients:
\[ c_n = \frac{2n+1}{2}\,\int_{-1}^i \mbox{sign}(x)\, P_n (x) \,{\text d}x = P_{n-1} (0) - P_{n+1} (0) , \qquad n=1,2,\ldots , \]
where \( P_{2k+1} (0) =0 \) and formulas for even indices was given previously. We now use Mathematica to calculate the coefficients up to arbitrary n, starting with n = 1.
Module[{i}, coef = {}; Do[coef = Append[coef, (LegendreP[i - 1, 0] - LegendreP[i + 1, 0])], {i, 1, 30}]]
The nth coefficient is then given by coef[[n]]. We sample a few values:
coef[[21]]
Out[34]= 180557/524288
coef[[22]]
Out[35]= 0
All the even coefficients are zero.

Now we define the mth partial sum, and then a graph of the mth partial sum, along with the original function. Because the coefficient of \( P_{0} (x)=1 \) is zero, we may start the sum over n with the n = 1 term. The given signum function is shown in blue and the partial sums in red.

legsum[m_, x_] := Module[{i}, Sum[N[coef[[i]]]*LegendreP[i, x], {i, 1, m}]]
legraph[m_] := Plot[{Sign[x], legsum[m, x]}, {x, -1.05, 1.05}, AxesLabel -> {"x", "Sign (x)"}, PlotRange -> {-1.5, 1.5}, PlotStyle -> {{RGBColor[0, 0, 1], Thickness[0.005]}, {RGBColor[1, 0, 0], Thickness [0.005]}}, PlotLabel -> Row[{"n=", PaddedForm[n, 3]}]]
legraph[30]
You may want to construct a sequence of partial sums which may be animated to see the convergence:
DynamicModule[{mangraph, i}, Do[mangraph[i] = legraph[i], {i, 1, 55, 2}]; Manipulate[mangraph[i], {i, 1, 51, 2}]]
You can use the slider to move through the graphs, and also to display a movie of the graphs. When you do this, you will see the familiar Gibbs overshoot/undershoot at the discontinuity. It is also observed a struggle to converge at the endpoints.
Module[{a}, coef = {}; Do[coef = Append[coef, (LegendreP[x - 1, 0] - LegendreP[x + 1, 0])], {x, 1, 100}]] legsum[m_, x_] := Module[{a}, -Sum[N[coef[[a]]]*LegendreP[a, x], {a, 1, m}]] leg[x_] = Piecewise[{{1, -1 <= x < 0}, {-1, 0 < x <= 1}}]; legraph[m_] := Plot[{leg[x], legsum[m, x]}, {x, -1, 1}, PlotRange -> {-1.5, 1.5}, PlotStyle -> {{RGBColor[0, 0, 1], Thickness[0.005]}, {RGBColor[1, 0, 0], Thickness[0.005]}}] legraph[10] legraph[100]

Example 3: The Dirac delta function admits the following Fourier--Legendre expansion:

\[ \delta (x-a) = \sum_{k\ge 0} \left( k + \frac{1}{2} \right) P_k (x)\, P_k (a) . \]
or
\[ \delta (x-a) = \frac{1}{2} \sum_{k\ge 0} \left( 2k+1 \right) P_k (x)\, P_k (a) . \]
So the partial sums
\[ \delta_N (x-a) = \frac{1}{2} \sum_{k= 0}^{N} \left( 2k+1 \right) P_k (x)\, P_k (a) \qquad \mbox{for} \quad |x| < 1, \]
have the property:
\[ \lim_{N\to \infty} \int_{-1}^1 {\text d}x \, \delta_N (x-a) \,f(x) = f(a) . \]
   ■

 

Return to Mathematica page
Return to the main page (APMA0340)
Return to the Part 1 Matrix Algebra
Return to the Part 2 Linear Systems of Ordinary Differential Equations
Return to the Part 3 Non-linear Systems of Ordinary Differential Equations
Return to the Part 4 Numerical Methods
Return to the Part 5 Fourier Series
Return to the Part 6 Partial Differential Equations
Return to the Part 7 Special Functions