---
title: "Sturm-Liouville theory"
firstLetter: "S"
publishDate: 2021-02-23
categories:
- Mathematics
- Physics

date: 2021-02-23T08:52:28+01:00
draft: false
markup: pandoc
---

# Sturm-Liouville theory

**Sturm-Liouville theory** defines the analogue of Hermitian matrix
eigenvalue problems for linear second-order ODEs.

It states that, given suitable boundary conditions, any linear
second-order ODE can be rewritten using the **Sturm-Liouville operator**,
and that the corresponding eigenvalue problem, known as a
**Sturm-Liouville problem**, will give real eigenvalues and a complete set
of eigenfunctions.


## General operator

Consider the most general form of a second-order linear
differential operator $\hat{L}$, where $p_0(x)$, $p_1(x)$, and $p_2(x)$
are real functions of $x \in [a,b]$ which are non-zero for all $x \in ]a, b[$:

$$\begin{aligned}
    \hat{L} \{u(x)\} = p_0(x) u''(x) + p_1(x) u'(x) + p_2(x) u(x)
\end{aligned}$$

We now define the **adjoint** or **Hermitian** operator
$\hat{L}^\dagger$ analogously to matrices:

$$\begin{aligned}
    \braket*{f}{\hat{L} g}
    = \braket*{\hat{L}^\dagger f}{g}
\end{aligned}$$

What is $\hat{L}^\dagger$, given the above definition of $\hat{L}$?
We start from the inner product $\braket*{f}{\hat{L} g}$:

$$\begin{aligned}
    \braket*{f}{\hat{L} g}
    &= \int_a^b f^*(x) \hat{L}\{g(x)\} \dd{x}
    = \int_a^b (f^* p_0) g'' + (f^* p_1) g' + (f^* p_2) g \dd{x}
    \\
    &= \big[ (f^* p_0) g' + (f^* p_1) g \big]_a^b - \int_a^b (f^* p_0)' g' + (f^* p_1)' g - (f^* p_2) g \dd{x}
    \\
    &= \big[ f^* \big( p_0 g' \!+\! p_1 g \big) \!-\! (f^* p_0)' g \big]_a^b + \int_a^b \! \big( (f p_0)'' - (f p_1)' + (f p_2) \big)^* g \dd{x}
    \\
    &= \big[ f^* \big( p_0 g' + (p_1 - p_0') g \big) - (f^*)' p_0 g \big]_a^b + \int_a^b \big( \hat{L}^\dagger\{f\} \big)^* g \dd{x}
\end{aligned}$$

We now have an expression for $\hat{L}^\dagger$, but are left with an
annoying boundary term:

$$\begin{aligned}
    \braket*{f}{\hat{L} g}
    &= \big[ f^* \big( p_0 g' + (p_1 - p_0') g \big) - (f^*)' p_0 g \big]_a^b + \braket*{\hat{L}^\dagger f}{g}
\end{aligned}$$

To fix this,
let us demand that $p_1(x) = p_0'(x)$ and that
$[p_0(f^* g' - (f^*)' g)]_a^b = 0$, leaving:

$$\begin{aligned}
    \braket*{f}{\hat{L} g}
    &= \big[ p_0 \big( f^* g' - (f^*)' g \big) \big]_a^b + \braket{\hat{L}^\dagger f}{g}
    = \braket*{\hat{L}^\dagger f}{g}
\end{aligned}$$

Using the aforementioned restriction $p_1(x) = p_0'(x)$,
we then take a look at the definition of $\hat{L}^\dagger$:

$$\begin{aligned}
    \hat{L}^\dagger \{f\}
    &= (p_0 f)'' - (p_1 f)' + (p_2 f)
    \\
    &= p_0 f'' + (2 p_0' - p_1) f' + (p_0'' - p_1' + p_2) f
    \\
    &= p_0 f'' + p_0' f' + p_2 f
    \\
    &= (p_0 f')' + p_2 f
\end{aligned}$$

The original operator $\hat{L}$ reduces to the same form,
so it is **self-adjoint**:

$$\begin{aligned}
    \hat{L} \{f\}
    &= p_0 f'' + p_0' f' + p_2 f
    = (p_0 f')' + p_2 f
    = \hat{L}^\dagger \{f\}
\end{aligned}$$

Consequently, every such second-order linear operator $\hat{L}$ is self-adjoint,
as long as it satisfies the constraints $p_1(x) = p_0'(x)$ and $[p_0 (f^* g' - (f^*)' g)]_a^b = 0$.

Let us ignore the latter constraint for now (it will return later),
and focus on the former: what if $\hat{L}$ does not satisfy $p_0' \neq p_1$?
We multiply it by an unknown $p(x) \neq 0$, and divide by $p_0(x) \neq 0$:

$$\begin{aligned}
    \frac{p(x)}{p_0(x)} \hat{L} \{u\} = p(x) u'' + p(x) \frac{p_1(x)}{p_0(x)} u' + p(x) \frac{p_2(x)}{p_0(x)} u
\end{aligned}$$

We now define $q(x)$,
and demand that the derivative $p'(x)$ of the unknown $p(x)$ satisfies:

$$\begin{aligned}
    q(x) = p(x) \frac{p_2(x)}{p_0(x)}
    \qquad
    p'(x) = p(x) \frac{p_1(x)}{p_0(x)}
\end{aligned}$$

The latter is a differential equation for $p(x)$, which we solve by integration:

$$\begin{gathered}
    \frac{p_1(x)}{p_0(x)} = \frac{1}{p(x)} \dv{p}{x}
    \quad \implies \quad
    \frac{p_1(x)}{p_0(x)} \dd{x} = \frac{1}{p(x)} \dd{p}
    \\
    \implies \quad
    \int_a^x \frac{p_1(\xi)}{p_0(\xi)} \dd{\xi} = \int_{p(a)}^{p(x)} \frac{1}{f} \dd{f}
    = \ln\Big( \frac{p(x)}{p(a)} \Big)
    \\
    \implies \quad
    p(x) = p(a) \exp\!\Big( \int_a^x \frac{p_1(\xi)}{p_0(\xi)} \dd{\xi} \Big)
\end{gathered}$$

Now that we have $p(x)$ and $q(x)$, we can define a new operator $\hat{L}_p$ as follows:

$$\begin{aligned}
    \hat{L}_p \{u\}
    = \frac{p}{p_0} \hat{L} \{u\}
    = p u'' + p' u' + q u
    = (p u')' + q u
\end{aligned}$$

This is the self-adjoint form from earlier!
So even if $p_0' \neq p_1$, any second-order linear operator with $p_0(x) \neq 0$
can easily be put in self-adjoint form.

This general form is known as the **Sturm-Liouville operator** $\hat{L}_{SL}$,
where $p(x)$ and $q(x)$ are non-zero real functions of the variable $x \in [a,b]$:

$$\begin{aligned}
    \boxed{
        \hat{L}_{SL} \{u(x)\}
        = \frac{d}{dx}\Big( p(x) \frac{du}{dx} \Big) + q(x) u(x)
        = \hat{L}_{SL}^\dagger \{u(x)\}
    }
\end{aligned}$$


## Eigenvalue problem

A **Sturm-Liouville problem** (SLP) is analogous to a matrix eigenvalue problem,
where $w(x)$ is a real weight function, $\lambda$ is the **eigenvalue**,
and $u(x)$ is the corresponding **eigenfunction**:

$$\begin{aligned}
    \boxed{
        \hat{L}_{SL}\{u(x)\} = - \lambda w(x) u(x)
    }
\end{aligned}$$

Necessarily, $w(x) > 0$ except in isolated points, where $w(x) = 0$ is allowed;
the point is that any inner product $\braket{f}{w g}$ may never be zero due to $w$'s fault.
Furthermore, the convention is that $u(x)$ cannot be trivially zero.

In our derivation of $\hat{L}_{SL}$,
we removed a boundary term to get self-adjointness.
Consequently, to have a valid SLP, the boundary conditions for
$u(x)$ must be as follows, otherwise the operator cannot be self-adjoint:

$$\begin{aligned}
    \Big[ p(x) \big( u^*(x) u'(x) - (u'(x))^* u(x) \big) \Big]_a^b = 0
\end{aligned}$$

There are many boundary conditions (BCs) which satisfy this requirement.
Some notable ones are listed here non-exhaustively:

+ **Dirichlet BCs**: $u(a) = u(b) = 0$
+ **Neumann BCs**: $u'(a) = u'(b) = 0$
+ **Robin BCs**: $\alpha_1 u(a) + \beta_1 u'(a) = \alpha_2 u(b) + \beta_2 u'(b) = 0$ with $\alpha_{1,2}, \beta_{1,2} \in \mathbb{R}$
+ **Periodic BCs**: $p(a) = p(b)$, $u(a) = u(b)$, and $u'(a) = u'(b)$
+ **Legendre "BCs"**: $p(a) = p(b) = 0$

Once this requirement is satisfied, Sturm-Liouville theory gives us
some very useful information about $\lambda$ and $u(x)$.
From the definition of an SLP, we know that, given two arbitrary (and possibly identical)
eigenfunctions $u_n$ and $u_m$, the following must be satisfied:

$$\begin{aligned}
    0 = \hat{L}_{SL}\{u_n\} + \lambda_n w u_n = \hat{L}_{SL}\{u_m^*\} + \lambda_m^* w u_m^*
\end{aligned}$$

We subtract these expressions, multiply by the eigenfunctions, and integrate:

$$\begin{aligned}
    0
    &= \int_a^b u_m^* \big(\hat{L}_{SL}\{u_n\} + \lambda_n w u_n\big) - u_n \big(\hat{L}_{SL}\{u_m^*\} + \lambda_m^* w u_m^*\big) \:dx
    \\
    &= \int_a^b u_m^* \hat{L}_{SL}\{u_n\} - u_n \hat{L}_{SL}\{u_m^*\} + u_n u_m^* w (\lambda_n - \lambda_m^*) \:dx
\end{aligned}$$

Rearranging this a bit reveals that these are in fact three inner products:

$$\begin{aligned}
    \int_a^b u_m^* \hat{L}_{SL}\{u_n\} - u_n \hat{L}_{SL}\{u_m^*\} \:dx
    &= (\lambda_m^* - \lambda_n) \int_a^b u_n u_m^* w \:dx
    \\
    \braket*{u_m}{\hat{L}_{SL} u_n} - \braket*{\hat{L}_{SL} u_m}{u_n}
    &= (\lambda_m^* - \lambda_n) \braket{u_m}{w u_n}
\end{aligned}$$

The operator $\hat{L}_{SL}$ is self-adjoint by definition,
so the left-hand side vanishes, leaving us with:

$$\begin{aligned}
    0
    &= (\lambda_m^* - \lambda_n) \braket{u_m}{w u_n}
\end{aligned}$$

When $m = n$, the inner product $\braket{u_n}{w u_n}$ is real and positive
(assuming $u_n$ is not trivially zero, in which case it would be disqualified anyway).
In this case we thus know that $\lambda_n^* = \lambda_n$,
i.e. the eigenvalue $\lambda_n$ is real for any $n$.

When $m \neq n$, then $\lambda_m^* - \lambda_n$ may or may not be zero,
depending on the degeneracy. If there is no degeneracy, we
see that $\braket{u_m}{w u_n} = 0$, i.e. the eigenfunctions are orthogonal.

In case of degeneracy, manual orthogonalization is needed, but as it turns out,
this is guaranteed to be doable, using e.g. the [Gram-Schmidt method](/know/concept/gram-schmidt-method/).

In conclusion, **a Sturm-Liouville problem has real eigenvalues $\lambda$,
and all the corresponding eigenfunctions $u(x)$ are mutually orthogonal**:

$$\begin{aligned}
    \boxed{
        \braket{u_m(x)}{w(x) u_n(x)}
        = \braket{u_n}{w u_n} \delta_{nm}
        = A_n \delta_{nm}
    }
\end{aligned}$$

When you're solving a differential eigenvalue problem,
knowing that all eigenvalues are real is a *huge* simplification,
so it is always worth checking whether you're dealing with an SLP.

Another useful fact of SLPs is that they always
have an infinite number of discrete eigenvalues.
Furthermore, the eigenvalues always ascend to $+\infty$;
in other words, there always exists a *lowest* eigenvalue $\lambda_0 > -\infty$,
known as the **ground state**.


## Completeness

Not only are the eigenfunctions $u_n(x)$ of an SLP orthogonal, they
also form a **complete basis**, meaning that any well-behaved function $f(x)$ can be
expanded as a **generalized Fourier series** with coefficients $a_n$:

$$\begin{aligned}
    \boxed{
        f(x)
        = \sum_{n = 0}^\infty a_n u_n(x)
        \quad \mathrm{for}\: x \in ]a, b[
    }
\end{aligned}$$

This series will converge significantly faster if $f(x)$
satisfies the same BCs as $u_n(x)$. In that case the
expansion will even be valid for the inclusive interval $x \in [a, b]$.

To find an expression for the coefficients $a_n$,
we multiply the above generalized Fourier series by $w(x) u_m^*(x)$ for an arbitrary $m$:

$$\begin{aligned}
    f(x) w(x) u_m^*(x)
    &= \sum_{n = 0}^\infty a_n u_n(x) w(x) u_m^*(x)
\end{aligned}$$

By integrating we get inner products on both the left and the right:

$$\begin{aligned}
    \int_a^b f(x) w(x) u_m^*(x) \dd{x}
    &= \int_a^b \Big(\sum_{n = 0}^\infty a_n u_n(x) w(x) u_m^*(x)\Big) \dd{x}
    \\
    \braket{u_m}{w f}
    &= \sum_{n = 0}^\infty a_n \braket{u_m}{w u_n}
\end{aligned}$$

Because the eigenfunctions of an SLP are mutually orthogonal,
the summation disappears:

$$\begin{aligned}
    \braket{u_m}{w f}
    &= \sum_{n = 0}^\infty a_n \braket{u_m}{w u_n}
    = \sum_{n = 0}^\infty a_n A_n \delta_{nm}
    = a_m A_m
\end{aligned}$$

After isolating this for $a_n$, we see that
the coefficients are given by the projection of the target
function $f(x)$ onto the normalized eigenfunctions $u_n(x) / A_n$:

$$\begin{aligned}
    \boxed{
        a_n
        = \frac{\braket{u_n}{w f}}{A_n}
        = \frac{\braket{u_n}{w f}}{\braket{u_n}{w u_n}}
    }
\end{aligned}$$

As a final remark, we can see something interesting
by rearranging the generalized Fourier series
after inserting the expression for $a_n$:

$$\begin{aligned}
    f(x)
    &= \sum_{n = 0}^\infty \frac{1}{A_n} \braket{u_n}{w f} u_n(x)
    = \int_a^b \Big(\sum_{n = 0}^\infty \frac{1}{A_n} u_n^*(\xi) w(\xi) f(\xi) u_n(x) \Big) \dd{\xi}
    \\
    &= \int_a^b f(\xi) \Big(\sum_{n = 0}^\infty \frac{1}{A_n} u_n^*(\xi) w(\xi) u_n(x) \Big) \dd{\xi}
    %= \int_a^b f(\xi) \delta(x - \xi) \dd{\xi}
\end{aligned}$$

Upon closer inspection, the parenthesized summation
must be the [Dirac delta function](/know/concept/dirac-delta-function/) $\delta(x)$
for the integral to work out.
This is in fact the underlying requirement for completeness:

$$\begin{aligned}
    \boxed{
        \sum_{n = 0}^\infty \frac{1}{A_n} u_n^*(\xi) w(\xi) u_n(x) = \delta(x - \xi)
    }
\end{aligned}$$



## References
1.  O. Bang,
    *Applied mathematics for physicists: lecture notes*, 2019,
    unpublished.