zenzicubi.co/posts/chebyshev/extra/index.qmd

---
title: "Generating Polynomials Extra: Legendary"
description: |
  "Another series of related, orthogonal polynomials"
format:
  html:
    html-math-method: katex
date: "2021-08-27"
date-modified: "2025-06-25"
categories:
  - algebra
  - generating functions
---

<style>
.figure-img {
  max-width: 512px;
  object-fit: contain;
  height: 100%;
}

.figure-img.wide {
  max-width: 768px;
}
</style>


In the [previous](../1) [two](../2) posts, I made clear the ties between the
  Chebyshev polynomials, constructibile polygons, and graph theory.
At this point, the benefits of using generating functions to reason about
  an infinite family of polynomials have hopefully been made manifest.

In this post, I will construct another family with different properties,
  but one that is still related to the Chebyshev polynomials.

Rather than taking the happy path, this post will focus less on taking the correct steps
  and more on 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 a list of coefficients with the following operations:

- Addition between two polynomials, implying
- Scalar multiplication of a polynomial
- Multiplication between two polynomials
- Evaluation at a number
    - Said "numbers", are the same type as the coefficients, e.g.,
      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
  we can define a dot product on them.
Again, from a data perspective, a dot product must be calculated from information
  in *both* polynomials and must have a value of the same type as the coefficients of the polynomial.

Naively, we could define one by using the already-defined multiplication on two polynomials,
  then evaluating the product at some value.
In actuality, the most typical choice is to *integrate* the product of two polynomials[^1] over a
  selected interval.
The simplest interval to consider is from -1 to 1, as this can be used to exploit symmetry for odd integrands.

[^1]: Sometimes a third term is included in the product, termed a *weight function*.
    The simplest weight function is 1, and is suitable for the family of polynomials we'll discuss shortly.

$$
P(z) \cdot Q(z) = \int_{-1}^1 P(z)Q(z)dz
$$

If the product of *P* and *Q* is an odd polynomial, they are trivially 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 ~ \left(\ = {2 \over 1}\right)
  \\
  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

- has degree 2,
- is orthogonal to both 1 and *z*,
- to continue the pattern, has norm 2/5?

Using a generic quadratic polynomial $p_2(z) = a z^2 + b z + c$, the second and third points
  give us constraints with which we can solve for *a*, *b*, and *c*.

$$
\begin{align*}
  \langle 1, p_2 \rangle
    &= \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
  \\
  &\implies a = -3c
  \\[14pt]
  \langle z, p_2 \rangle
    &= \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
  \\
  &\implies b = 0
  \\[14pt]
  \langle p_2, p_2 \rangle
    &= \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 we can see that the new polynomial is completely even, like $1$ which preceded it.
The final integral, though not challenging, requires some additional steps to resemble an equation
  in *a* and *c* to solves the system.

$$
\begin{align*}
  \langle p_2, p_2 \rangle
    &= \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
  \\
  &= {18 c^2 \over 5} - {12 c^2 \over 3} + 2c^2
  \\
  &= {18 \cdot 3 - 12 \cdot 5 + 2 \cdot 15 \over 15} c^2
    = {24 \over 15} c^2 = {8 \over 5} c^2
    = {2 \over 5}
  \\
  &\implies 4c^2 = 1 \qquad c = \pm {1 \over 2}
  \\
  &\implies 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 instead used constraint requiring less algebra, such as by forcing
  $p_2(1) = 1$ (though this turns out to be equivalent to the one already given).

Since we now have three polynomials, we can now build a new cubic polynomial
  with the same properties in a similar manner.
Doing so ad infinitum would build an entire family of mutually orthogonal polynomials.
This family is called the
  [Legendre polynomials](https://en.wikipedia.org/wiki/Legendre_polynomials),
  named after Adrien-Marie Legendre.
Intriguingly, only a few surviving portraits of this mathematician are known, one of which is a fairly humorous caricature.

<!-- TODO: better wikimedia imports -->
<div class="quarto-figure quarto-figure-center">
  <figure class="figure">
  <a title="Julien-Léopold Boilly, Public domain, via Wikimedia Commons" href="https://commons.wikimedia.org/wiki/File:Legendre.jpg">
    <img class="img-fluid figure-img" alt="Legendre" src="https://upload.wikimedia.org/wikipedia/commons/0/03/Legendre.jpg?20210711144200">
  </a>

  <figcaption>
    Public domain image retrieved from Wikimedia
  </figcaption>
  </figure>
</div>

An orthogonal family of polynomials is useful when doing numeric computations involving integrals.
Expressing a (scaled and shifted) general polynomial in an orthogonal basis not only
  simplifies the math, but can prevent errors from accumulating.
These polynomials also have applications in electrical engineering.

However, calculating them naively will require more free coefficients, more integrals,
  and most importantly, more effort.
Is there a way to bypass so much math? How did Legendre do it?


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

Of course there's an easier way.
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*

Instead of dealing with the polynomials individually, we can, in the same manner
  as with the Chebyshev polynomials, use a generating function to gather all of them into a single expression.
As a reminder (and for the sake of consistency with the previous post),
  the polynomial variable will be *z*, and the series variable will be *x*.

We don't know what the generating function 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 conditions, which exist for the integral
  of 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_{\substack{
      n = 0 \\
      r + s = n \\
      r \neq s
    }}^\infty p_r p_s x^n
    + \sum_{n = 0}^\infty p_n^2 x^{2n}
\end{align*}
$$

This expression can be integrated over *z* on the boundary defined by the inner product.
By linearity, 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) dz
  \\
  &= \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} + 0
\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*}
$$


### A wrong guess

Since we integrated over *z*, it does not appear in the resulting expression.
On the other hand, as a difference of two terms, the expression resembles an evaluated integral in *z*,
  so it's possible to conjecture what the antiderivative looks like:

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

Now we can take the derivative of this expression with respect to *z* to get
  a possible value for $F^2$, and hence *F*.

$$
\begin{align*}
  \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*}
$$

This expression is clearly wrong.
It can be interpreted as a power series in *zx*, meaning that as a series in *x*, there is
  only one *z* in a coefficient, with power equal to the power of *x*.
For example, 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$.


### The correction

Using logarithm identities, we can take the term inside to a higher power and divide out by it.
Doing so, then differentiating as before, yields a slightly different expression.

$$
\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*}
$$

We can justify this manipulation by noticing that we start with two polynomials (1 and *z*),
  so the denominator must be quadratic in *x* to fit with the two initial values.

We can then examine first two values of the series by evaluating the series and its derivative at 0.

$$
\begin{align*}
  f_0(z) &= F(0; z)
  \\
  &= {1 \over \sqrt{1 + 2z \cdot 0 + 0^2} } = 1
  \\[12pt]
  f_1(z)
    &= \partial_x F(x; z) |_{x = 0}
  \\
  &= \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*}
$$

While the 0th term agrees with $p_0$, the first term has the opposite sign of $p_1$.
We can amend our expression for *F* by replacing *z* with -*z*.
Then, we can integrate it with respect to *z* to verify it has the same properties.

$$
\begin{align*}
  F(x; z)
    &\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*}
$$

This is the exact same expression we obtained from 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 expression for *F* is correct, and the generating function
  for the Legendre polynomials is indeed

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

Doesn't this look somewhat familiar?
It is simply the square root of the expression which generates the Chebyshev polynomials
  of the second kind (with zeroth term 1).

$$
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}
    = \partial_x {1 \over \sqrt{1 - 2zx + x^2} }
  \\
  &= {z - x \over ( 1 - 2zx + x^2 )^{3/2} }
    = {z - x \over 1 - 2zx + x^2 } F(x; z)
  \\
  (z - x) F(x; z) &= (1 - 2zx + x^2) \partial_x F(x; z)
  \\
  &\implies (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 easiest thing to do is distribute each product and adjust the bounds of the summation to start at 1.
It's also convenient to sum over $x^n$ rather than other powers.

$$
\begin{gather*}
  \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.
\end{gather*}
$$

Now the only problem which remains are the elements which are *not* in a sum, $zp_0$ and $p_1$.
But these are both equal to *z*, and they are on opposite sides of the equation, so they cancel out.
With the sums aligned, they can all be dropped from the expression, giving an explicit recurrence.

$$
\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}u
  \\
  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,
  Legendre 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.