---
format:
  html:
    html-math-method: katex
---


Generating Polynomials Extra: Legendary
=======================================

In the [previous]() [two]() posts, I made clear the ties between the Chebyshev polynomials and polygonal constructibility and graph theory, the benefits of having an infinite family of polynomials are manifold. In this post, I will construct another family with different properties, but one that is still related to the Chebyshev polynomials.

The format of this post will be slightly different: it will have a less concrete goal and be geared more toward following mathematical intuition.


Orthogonality
-------------

Another very basic question we can ask is: what exactly is a polynomial? Formally (or by viewing it as a container for data), it is no more than a list of numbers with the following operations:

- Addition between two polynomials, implying
- Scalar multiplication of a polynomial
- Multiplication between two polynomials
- Evaluation at a number

By "number", I mean an object of the same type as the coefficients, which can be a natural, an integer, a rational, etc. The first two properties above and existence of the zero polynomial imply that polynomials can be considered a [vector space](https://en.wikipedia.org/wiki/Vector_space).

Another thing they have in common with typical Cartesian vectors is that a dot product can be defined on them. Again from a data perspective, a dot product must be calculated from information in *both* polynomials and have a value whose type is the same as the coefficients of the polynomial. Naively, this could be done by using the already-defined multiplication on two polynomials, then evaluating the product at some value.

In actuality, this is done by *integrating* the product of two polynomials. Sometimes a third expression is included in the product, termed a *weight function*, but the simplest weight function is 1. Similarly, the simplest interval is from -1 to 1, as this can be used to exploit symmetry for odd integrands.

$$
\begin{gather}
P(z) \cdot_w Q(z) = \int_a^b P(z)Q(z)w(z)dz \\
P(z) \cdot Q(z) = \int_{-1}^1 P(z)Q(z)dz
\end{gather}
$$

If the product of *P* and *Q* is an odd polynomial, then their dot product is trivially 0, and they are orthogonal. For example, if $P(z) = 1$ and $Q(z) = z$, then the integrand is simply *z*; this polynomial is odd, so $1$ and $z$ are orthogonal. The inner products of $1$ and $z$ with themselves (or norms, by analogy with typical vectors) are:

$$
\begin{gather}
1 \cdot 1 = \int_{-1}^1 1^2 dz = \left.{z^{\vphantom{0}}}\right]_{-1}^1 = 1 -\ (-1) = 2 \\
z \cdot z = \int_{-1}^1 z^2 dz = \left.{z^3 \over 3}\right]_{-1}^1 = {1 \over 3} -\ \left(-{1 \over 3}\right) = {2 \over 3}
\end{gather}
$$

What if we were to add a third polynomial into the mix, which is orthogonal to both $1$ and $z$, has degree 2, and whose norm is, say, 2/5? To start, just use a generic quadratic polynomial $p_2(z) = a z^2 + b z + c$. Then enforce all three constraints:

$$
\begin{align*}
1 \cdot p_2 &= \int_{-1}^1 (a z^2 + b z + c) dz
= \left[ {a z^3 \over 3} + c z \right]_{-1}^1
= {2 \over 3}a + 2c = 0 \\
z \cdot p_2 &= \int_{-1}^1 (a z^3 + b z^2 + c z) dz
= \left[ {b z^3 \over 3} \right]_{-1}^1
= {2 \over 3}b = 0 \\
p_2 \cdot p_2 &= \int_{-1}^1 (a z^2 + c)^2 dz
= \int_{-1}^1 (a^2 z^4 + 2ac z^2 + c^2) dz = {2 \over 5}
\end{align*}
$$

Already it can be seen that *b* is 0, which implies the new polynomial is completely even, like $1$ which preceded it. The final integral, though not challenging, requires a fair number of steps to resemble a constraint on *a* and *c* which solves the system.

$$
\begin{align*}
\int_{-1}^1 (a^2 z^4 + 2ac z^2 + c^2) dz
=&~ \left[ a^2 {z^5 \over 5} + ( 2ac ) {z^3 \over 3} + c^2 z \right]_{-1}^1 \\
=&~ {2a^2 \over 5} + {4ac \over 3} + 2c^2
= {2a^2 \over 5} + {2ac \over 3} + {2ac \over 3} + 2c^2 \\
=&~ {2 \over 15} ( 3a^2 + 5ac ) + c \left( {2a \over 3} + 2c \right) \\
=&~ {2 \over 15} ( 3(-3c)^2 + 5(-3c)c ) + c (0)
= {2 \over 15}(27c^2 -\ 15c^2) \\
=&~ {2 \over 5} (4c^2) = {2 \over 5} \implies c = \pm {1 \over 2} \\
\implies& {2a \over 3} \pm 1 = 0 \qquad a = \mp {3 \over 2}
\end{align*}
$$

Preferring the leading coefficient a to be positive, this yields the polynomial $p_2(z) = {1 \over 2}(3z^2 -\ 1)$. We might have sought a different constraint which requires less algebra to evaluate, for example by forcing $p_2(1) = 1$ (though this turns out to be equivalent to the one already given).

This family of polynomials is called the [Legendre polynomials](https://en.wikipedia.org/wiki/Legendre_polynomials), named after Adrien-Marie Legendre. I bring up the man behind them in particular because only one surviving portrait of him exists, and it is a fairly humorous caricature.

::: {}
![]()
Public domain image retrieved from Wikimedia
:::

An entire family of mutually orthogonal polynomials has its uses. When doing numeric computation involving integrals, expressing a (scaled and shifted) general polynomial in the orthogonal basis not only simplifies the math, but can prevent errors from accumulating. They also have applications in electrical engineering.

Even with the alternative, it becomes clear that extending the family further will require more free coefficients, more integrals, and an overall bigger system. How did Legendre do it? Is there a way to bypass so much math?


Presupposition and Trial and Error
----------------------------------

Of course there is. What do we know about this problem? We're dealing with a sequence of polynomials, where

- A polynomial's degree uniquely identifies it in the sequence,
- Any two polynomials which are not the same are orthogonal,
- Any polynomial's norm is $2/(2n + 1)$, where *n* is the degree of the polynomial
- The first two terms are $1$ and $z$

In the same manner as with the Chebyshev polynomials, considering a sequence of polynomials is easiest if given as an expression whose Taylor series coefficients are the requisite polynomials. As a reminder, and for the sake of consistency with the previous post, the polynomial variable has been *z*, and the series variable will be *x*.

We don't know what this expression is, and we only know a few terms of the series

$$
\begin{align*}
F(x; z) &= \sum_{n=0}^\infty p_n(z) x^n \\
&= p_0(z) + p_1(z)x + p_2(z)x^2 + ... \\
&= 1 + zx + {1 \over 2}(3z^2 -\ 1) x^2 + ...
\end{align*}
$$

We want to exploit the norm and orthogonality condition, which are obtained by integrating the product of any two coefficients. As a single object *F*, it is very easy to get every possible product of two terms -- it is simply given by $F^2$. Despite the apparent difficulty in distributing over an infinite number of terms, the result can be expressed somewhat simply.

$$
\begin{align*}
F(x; z)^2 &= \left( \sum_{n=0}^\infty p_n x^n \right)^2
= (p_0 + p_1 x + p_2 x^2 + ...)(p_0 + p_1 x + p_2 x^2 + ...) \\
&= p_0^2 + (2p_0 p_1) x + (2p_0 p_2 + p_1^2) x^2 + (2p_0 p_3 + 2p_1 p_2) x^3 + ... \\
&= \sum_{\substack{n = 0 \\ r + s = n}}^\infty p_r p_s x^n
= \sum_{n = 0}^\infty p_n^2 x^{2n}
+ \sum_{\substack{n = 0 \\ r + s = n \\ r \neq s}}^\infty p_r p_s x^n
\end{align*}
$$

This expression can be integrated over *z* on the boundary defined by the inner product. The integral of the sum is the sum of the integrals since x as a series variable is totally independent of *z* (as well as [a separate kind of object](https://en.wikipedia.org/wiki/Indeterminate_%28variable%29)).

$$
\begin{align*}
\int_{-1}^1 F(x; z)^2 dz
&= \int_{-1}^1 \left( \sum_{n = 0}^\infty p_n^2 x^{2n} \right) dz
+ \int_{-1}^1 \left( \sum_{
\substack{n = 0 \\ r + s = n \\ r \neq s}}^\infty p_r p_s x^n \right) \\
&= \sum_{n = 0}^\infty \left( \int_{-1}^1 p_n^2 dz \right) x^{2n}
+ \sum_{\substack{n = 0 \\ r + s = n \\ r \neq s}}^\infty
\left( \int_{-1}^1 p_r p_s dz \right) x^n \\
&= \sum_{n = 0}^\infty {2 \over 2n + 1} x^{2n}
\end{align*}
$$

When multiplied by *x*, this sum very much looks like the integral of another expression,

$$
\begin{align*}
x \int_{-1}^1 F(x; z)^2 dz
&= \sum_{n = 0}^\infty {2 \over 2n + 1} x^{2n+1}
= \int \left( \sum_{n = 0}^\infty 2 x^{2n} \right) dx \\
&= \int {2dx \over 1 -\ x^2} = \int { \left( {A \over 1 + x} + {B \over 1 -\ x} \right) } \\
&= A\ln(1+x) -\ B\ln(1 -\ x) \\
2 &= A(1 -\ x) + B(1 + x)
\implies \stackrel{x = -1}{A = 1}, \quad \stackrel{x = 1}{B = 1} \\
\end{align*}
$$

As a reminder, the goal of these manipulations is to find a closed expression which resembles the evaluation of integral on the LHS with respect to *z*. Fortunately, the partial fraction decomposition is not only simple, but resembles an evaluated integral.

Before continuing, it is imperative that *x* be divided out of both sides. This implies that the coefficient of *z* in the expression of the logarithm is a multiple of *x*. For example:

$$
\begin{align*}
\int_{-1}^1 F^2 dz
&= {1 \over x}( \ln(1+x) -\ \ln(1 -\ x) ) \\
&= {1 \over x} \left[ \ln(1 +\ zx)^{\vphantom{0}} \right]_{z = -1}^1 \\
\partial_z \ln(1 + zx) &= {x \over 1 + zx} \\
\implies F^2 &\stackrel{?}{=} {1 \over 1 + zx}
\quad F \stackrel{?}{=} {1 \over \sqrt{1 + zx} }
\end{align*}
$$

But this expression is clearly wrong. It can be interpreted as a power series in *zx*, rather than just in *x*. This series can also come from a binomial expansion in $1 + zx$, and its terms compared with the ones we already know. The second Legendre polynomial is ${1 \over 2}(3z^2 -\ 1)$, but it is impossible that this expression will be generated, since it contains a term other than a multiple of $z^2$.

Using logarithm identities, an alternate expression can be obtained

$$
\begin{align*}
\int_{-1}^1 F^2 dz
&= {1 \over x}( \ln(1+x) -\ \ln(1 -\ x) ) \\
&= {1 \over 2x}( \ln((1+x)^2) -\ \ln((1 -\ x)^2) )\\
&= {1 \over 2x}( \ln(1+ 2x + x^2) -\ {1 \over 2}\ln(1 -\ 2x + x^2) ) \\
&= {1 \over 2x} \left[ \ln(1 + 2zx + x^2) \right]_{z = -1}^1 \\
\partial_z \ln(1 + 2zx + x^2) &= {2x \over 1 + 2zx + x^2} \\
\implies F^2 &\stackrel{?}{=} {1 \over 1 + 2zx + x^2}
\quad F \stackrel{?}{=} {1 \over \sqrt{1 + 2zx + x^2} }
\end{align*}
$$

The justification for this move is that we start with two polynomials, $1$ and $z$, so the expression must be quadratic in *x* to fit with the two initial values. Examining the first term of the series (i.e., the coefficient of *x*) by taking its derivative, this expression is only *almost* correct.

$$
\begin{align*}
\left. \left( \partial_x (1 + 2zx + x^2)^{-1/2} \right) \right|_{x = 0} &=
\left. -{1 \over 2} (1 + 2zx + x^2)^{-3/2} (2z + 2x) \right|_{x = 0} \\
&= (1 + 0 + 0)^{-3/2} (-z -\ 0) \\
&= -z
\end{align*}
$$

This has the opposite sign of $z$, the first term in the series. However, the series does correctly generate the zeroth term 1. Thus, by replacing *z* with -*z* in the expression of *F*, one obtains a new expression, which can be integrated to verify its validity.

$$
\begin{align*}
F &\stackrel{?}{=} {1 \over \sqrt{1 -\ 2zx + x^2} } \\
\int_{-1}^1 {dz \over 1 -\ 2zx + x^2}
&= -{1 \over 2x}\left[ \ln( 1 -\ 2zx + x^2 ) \right]_{z=-1}^{1} \\
&= -{1 \over 2x}( \ln( 1 -\ 2x + x^2 ) -\ \ln( 1 + 2x + x^2 ) ) \\
&= -{1 \over 2x}( \ln( (1 -\ x)^2 ) -\ \ln( (1 + x)^2 ) ) \\
&= {1 \over x}( \ln( 1 + x ) -\ \ln( 1 -\ x ) )
\end{align*}
$$

Which is the exact same as the previous constraint given by the orthogonality condition. These previous hurdles are in part due to the non-uniqueness of the square root, but also because reversing integration over symmetric bounds is challenging.

This final expression is correct, and the generating function for the Legendre polynomials is

$$
F(x; z) = {1 \over \sqrt{1 -\ 2zx + x^2} }
$$

Isn't this form somewhat familiar? It is simply the square root of the expression which generates the Chebyshev polynomials of the second kind (without an offset).

$$
F(x; z) = \sqrt{B(x; z) \over x}
$$

The generating function approach is (at least according to Wikipedia) how Legendre approached this problem. From the math above, it certainly has a degree of elegance to it, albeit one which requires some apparent (but justified) leaps in logic.


Recurrence
----------

With a closed expression finally known, is it possible to obtain a recurrence? It is, but doing so is no easier than the previous steps. The result and process for doing so are outlined on Wikipedia, but I will do so explicitly here. Begin by differentiating *F* with respect to *x*, which shifts the series and multiplies each polynomial by its place in the sequence (i.e., its degree)

$$
\begin{align*}
F(x; z) &= \sum_{n=0}^\infty p_n(z)x^n
= {1 \over \sqrt{1 -\ 2zx + x^2} } \\
\partial_x F(x; z) &= \sum_{n=1}^\infty n p_n(z)x^{n-1}
= {z -\ x \over ( 1 -\ 2zx + x^2 )^{3/2} }
= {z -\ x \over 1 -\ 2zx + x^2 } F(x; z) \\
(z -\ x) \sum_{n=0}^\infty p_n x^n
&= (1 -\ 2zx + x^2) \sum_{n=1}^\infty n p_n x^{n-1}
\end{align*}
$$

With the sums re-exposed, the most difficult part is aligning like terms. The most natural range of summation starts at 1 (like the RHS), but the power of *x* is simply $x^n$ (like the LHS). Therefore, distributing over each sum and realigning indices yields:

$$
\text{LHS} \left \{ \begin{align*}
z\sum_{n=0}^\infty p_n x^n
&= z p_0 + z \sum_{n=1}^\infty p_n x^n
\\
x\sum_{n=0}^\infty p_n x^n
&= \sum_{n=0}^\infty p_n x^{n+1}
= \sum_{n=1}^\infty p_{n-1} x^n
\end{align*} \right .
\\
\text{RHS} \left \{ \begin{align*}
\sum_{n=1}^\infty n p_n x^{n-1}
&= \sum_{n=0}^\infty (n + 1) p_{n + 1} x^n
= p_1 + \sum_{n=1}^\infty (n + 1) p_{n + 1} x^n
\\
2zx\sum_{n=1}^\infty n p_n x^{n-1}
&= 2z\sum_{n=1}^\infty n p_n x^n \\
x^2\sum_{n=1}^\infty n p_n x^{n-1}
&= \sum_{n=1}^\infty n p_n x^{n+1}
= \sum_{n=0}^\infty n p_n x^{n+1}
= \sum_{n=1}^\infty (n -\ 1) p_{n-1} x^n
\end{align*}\right.
$$

Now the only problem which remains are the elements which are *not* in a sum, $zp_0$ and $p_1$. But these are both just *z*, and they are on opposite sides of the equation, so they cancel out. The remaining sums can be dropped entirely, and an explicit recurrence is formed.

$$
\begin{align*}
(z -\ x) \sum_{n=0}^\infty p_n x^n
&= (1 -\ 2zx + x^2) \sum_{n=1}^\infty n p_n x^{n-1} \\
z \sum_{n=1}^\infty p_n x^n -\ \sum_{n=1}^\infty p_{n-1} x^n
&= \sum_{n=1}^\infty (n + 1) p_{n + 1} x^n -\ 2z\sum_{n=1}^\infty n p_n x^n
+ \sum_{n=1}^\infty (n -\ 1) p_{n-1} x^n \\
z p_n -\ p_{n-1} &=
(n + 1) p_{n + 1} -\ 2z n p_n + (n -\ 1) p_{n-1} \\
-(n + 1) p_{n + 1} &=
-z p_n (n + 1) -\ 2z n p_n + (n -\ 1) p_{n-1} + p_{n-1} \\
(n + 1) p_{n + 1} &= (2n + 1)z p_n -\ n p_{n-1}
\end{align*}
$$

This form is called *Bonnet’s recursion formula*. Trying out this formula with $1$ and $z$ indeed produces the correct next term.

$$
\begin{align*}
(n + 1) p_{n + 1} &= (2n + 1)z p_n -\ n p_{n-1} \\
(1 + 1) p_{1 + 1} &= (2(1) + 1)z p_1 -\ 1 p_{1-1} \\
2 p_2 &= 3z(z) -\ p_0 \\
p_2 &= {1 \over 2}(3z^2 -\ 1)
\end{align*}
$$

The recursion formula also implies something interesting: like the Chebyshev polynomials of the second kind, all polynomials are either totally even or totally odd.

***

With the Chebyshev polynomials, the goal was to manipulate a recurrence relation into a generating function. However, the Legendre polynomials require the opposite; one first starts by assuming they have the generating function. Only through clever manipulations and experimentation can an expression be recovered, and from the expression, a recurrence.
