Sturm-Liouville Theory

Sturm-Liouville theory provides the rigorous mathematical framework underlying eigenfunction expansions in physics. It unifies the theory of Fourier series, Bessel functions, Legendre polynomials, and quantum mechanical wave functions into a single elegant structure: self-adjoint differential operators with real eigenvalues and complete orthogonal eigenfunctions.

1. Self-Adjoint Operators

The Sturm-Liouville operator is a second-order linear differential operator of the form:

$$\mathcal{L}[y] = \frac{1}{w(x)}\left[-\frac{d}{dx}\!\left(p(x)\frac{dy}{dx}\right) + q(x)\,y\right]$$

where $p(x) > 0$, $w(x) > 0$ (the weight function), and $q(x)$ is real on the interval $[a, b]$. The Sturm-Liouville eigenvalue problem is:

$$-\frac{d}{dx}\!\left(p(x)\frac{dy}{dx}\right) + q(x)\,y = \lambda\, w(x)\, y$$

subject to boundary conditions at $x = a$ and $x = b$.

1.1 The Self-Adjointness Property

An operator $\mathcal{L}$ is self-adjoint (Hermitian) with respect to the weighted inner product if:

$$\langle u, \mathcal{L}v \rangle_w = \langle \mathcal{L}u, v \rangle_w$$

where the weighted inner product is $\langle u, v \rangle_w = \int_a^b u^*(x)\,v(x)\,w(x)\,dx$.

Proof of self-adjointness: Consider $\langle u, \mathcal{L}v \rangle_w - \langle \mathcal{L}u, v \rangle_w$. Using the operator definition (without the $1/w$ factor):

$$\int_a^b u^*\!\left[-(pv')' + qv\right]dx - \int_a^b \left[-(pu^*)' + qu^*\right]v\,dx = -\int_a^b \left[u^*(pv')' - (pu^*)'v\right]dx$$

Integrating by parts (Green's identity):

$$= -\left[p(x)\!\left(u^*v' - u^{*\prime}v\right)\right]_a^b$$

This boundary term vanishes for the standard boundary conditions:

Dirichlet: $y(a) = y(b) = 0$
Neumann: $y'(a) = y'(b) = 0$
Periodic: $y(a) = y(b)$, $y'(a) = y'(b)$
Mixed (Robin): $\alpha y(a) + \beta y'(a) = 0$
Singular: $p(a) = 0$ or $p(b) = 0$ (e.g., Legendre equation at $x = \pm 1$)

2. Eigenvalue Problems: Key Theorems

Self-adjointness leads to remarkable properties of the eigenvalue problem$\mathcal{L}y_n = \lambda_n y_n$:

2.1 Real Eigenvalues

Theorem: All eigenvalues $\lambda_n$ are real.

Proof: Suppose $\mathcal{L}y = \lambda y$. Taking the inner product with $y$:

$$\langle y, \mathcal{L}y \rangle_w = \lambda \langle y, y \rangle_w$$

Taking the complex conjugate and using self-adjointness:

$$\overline{\langle y, \mathcal{L}y \rangle_w} = \langle \mathcal{L}y, y \rangle_w = \langle y, \mathcal{L}y \rangle_w = \bar{\lambda}\langle y, y \rangle_w$$

Therefore $(\lambda - \bar{\lambda})\langle y, y \rangle_w = 0$. Since$\langle y, y \rangle_w > 0$ for non-trivial $y$, we conclude$\lambda = \bar{\lambda}$, i.e., $\lambda \in \mathbb{R}$.

2.2 Orthogonality of Eigenfunctions

Theorem: Eigenfunctions corresponding to distinct eigenvalues are orthogonal with respect to the weight function.

Proof: Let $\mathcal{L}y_n = \lambda_n y_n$ and$\mathcal{L}y_m = \lambda_m y_m$ with $\lambda_n \neq \lambda_m$. Then:

$$\langle y_m, \mathcal{L}y_n \rangle_w = \lambda_n \langle y_m, y_n \rangle_w$$

$$\langle \mathcal{L}y_m, y_n \rangle_w = \lambda_m \langle y_m, y_n \rangle_w$$

By self-adjointness, the left sides are equal, so$(\lambda_n - \lambda_m)\langle y_m, y_n \rangle_w = 0$. Since$\lambda_n \neq \lambda_m$:

$$\langle y_m, y_n \rangle_w = \int_a^b y_m^*(x)\, y_n(x)\, w(x)\,dx = 0$$

2.3 Discrete Spectrum and Ordering

For a regular Sturm-Liouville problem on a finite interval, the eigenvalues form a countably infinite increasing sequence:

$$\lambda_1 < \lambda_2 < \lambda_3 < \cdots \to \infty$$

The $n$-th eigenfunction $y_n(x)$ has exactly $n-1$ zeros in the open interval $(a, b)$. This oscillation theorem is analogous to the quantum mechanical result that excited states have more nodes.

3. Completeness of Eigenfunctions

The eigenfunctions $\{y_n\}$ of a Sturm-Liouville problem form a completeset in $L^2_w[a,b]$: any square-integrable function can be expanded as:

$$f(x) = \sum_{n=1}^{\infty} c_n\, y_n(x), \quad c_n = \frac{\langle y_n, f \rangle_w}{\langle y_n, y_n \rangle_w} = \frac{\int_a^b y_n^*(x)\,f(x)\,w(x)\,dx}{\int_a^b |y_n(x)|^2\,w(x)\,dx}$$

Convergence is in the $L^2_w$ sense (mean-square), meaning:

$$\lim_{N\to\infty} \int_a^b \left|f(x) - \sum_{n=1}^{N} c_n\, y_n(x)\right|^2 w(x)\,dx = 0$$

3.1 The Completeness Relation

Completeness can be expressed as a resolution of the identity. For normalized eigenfunctions ($\langle y_n, y_n \rangle_w = 1$):

$$\sum_{n=1}^{\infty} y_n^*(x')\, y_n(x) = \frac{\delta(x - x')}{w(x)}$$

This is the closure relation, analogous to $\sum_n |n\rangle\langle n| = \hat{I}$ in Dirac notation. It guarantees that the eigenfunction basis spans the entire function space.

3.2 Generalized Fourier Series

Every eigenfunction expansion in physics is a special case of the Sturm-Liouville completeness theorem:

Problem	$p(x)$	$w(x)$	Eigenfunctions
Fourier series	$1$	$1$	$\sin(n\pi x/L)$, $\cos(n\pi x/L)$
Legendre	$1-x^2$	$1$	$P_\ell(x)$
Bessel	$x$	$x$	$J_m(\alpha_{mn}x/a)$
Hermite (QM)	$e^{-x^2}$	$e^{-x^2}$	$H_n(x)e^{-x^2/2}$
Laguerre (QM)	$xe^{-x}$	$e^{-x}$	$L_n^\alpha(x)$

4. Green's Function Expansion

The Green's function for the Sturm-Liouville operator satisfies:

$$\mathcal{L}G(x, x') = \frac{\delta(x - x')}{w(x)}$$

with homogeneous boundary conditions. Using the completeness of eigenfunctions, we can expand both $G$ and $\delta/w$:

4.1 Eigenfunction Expansion of the Green's Function

Derivation: Expand $G(x, x') = \sum_n a_n(x')\, y_n(x)$. Applying$\mathcal{L}$:

$$\sum_n a_n(x')\,\lambda_n\, y_n(x) = \frac{\delta(x - x')}{w(x)} = \sum_n y_n^*(x')\, y_n(x)$$

where we used the completeness relation in the last step (for normalized eigenfunctions). Comparing coefficients: $a_n(x') = y_n^*(x') / \lambda_n$. Therefore:

$$G(x, x') = \sum_{n=1}^{\infty} \frac{y_n(x)\, y_n^*(x')}{\lambda_n}$$

This is known as the bilinear expansion (or Mercer expansion) of the Green's function. It converges provided $\lambda_n \neq 0$ for all $n$.

4.2 Solving Inhomogeneous Equations

The solution to the inhomogeneous equation $\mathcal{L}u = f$ is:

$$u(x) = \int_a^b G(x, x')\, f(x')\, w(x')\,dx' = \sum_{n=1}^{\infty} \frac{f_n}{\lambda_n}\, y_n(x)$$

where $f_n = \langle y_n, f \rangle_w$. This is the eigenfunction method for solving boundary value problems — exactly analogous to diagonalizing a matrix.

4.3 Modified Green's Function

If $\lambda = 0$ is an eigenvalue (e.g., the Neumann problem for Laplace's equation), the Green's function does not exist in the usual sense. The modified Green's functionexcludes the zero-mode:

$$G_m(x, x') = \sum_{n=2}^{\infty} \frac{y_n(x)\, y_n^*(x')}{\lambda_n}$$

A solution then exists only if $f$ is orthogonal to the zero-mode eigenfunction (the Fredholm alternative).

5. Connection to Quantum Mechanics

The time-independent Schrodinger equation is a Sturm-Liouville problem:

$$-\frac{\hbar^2}{2m}\frac{d^2\psi}{dx^2} + V(x)\psi = E\psi$$

This is precisely the Sturm-Liouville form with $p(x) = \hbar^2/(2m)$,$q(x) = V(x)$, $w(x) = 1$, and $\lambda = E$.

5.1 The Quantum Postulates as SL Theory

The fundamental postulates of quantum mechanics are direct consequences of Sturm-Liouville theory:

Real eigenvalues $\to$ Measured energies are real
Orthogonality $\to$ Distinct energy eigenstates are orthogonal: $\langle \psi_n | \psi_m \rangle = \delta_{nm}$
Completeness $\to$ Any state can be expanded: $|\Psi\rangle = \sum_n c_n |\psi_n\rangle$
Oscillation theorem $\to$ Ground state has no nodes, $n$-th excited state has $n$ nodes

5.2 The Harmonic Oscillator

The quantum harmonic oscillator $V(x) = \frac{1}{2}m\omega^2 x^2$ gives the Hermite equation in Sturm-Liouville form. The eigenfunctions are:

$$\psi_n(x) = \left(\frac{m\omega}{\pi\hbar}\right)^{1/4} \frac{1}{\sqrt{2^n n!}} H_n\!\left(\sqrt{\frac{m\omega}{\hbar}}x\right) e^{-m\omega x^2/(2\hbar)}$$

with eigenvalues $E_n = \hbar\omega(n + 1/2)$. The orthogonality relation$\langle\psi_n|\psi_m\rangle = \delta_{nm}$ is precisely the Sturm-Liouville orthogonality with weight $w(x) = 1$.

5.3 The Hydrogen Atom

The radial Schrodinger equation for hydrogen, with $V(r) = -e^2/(4\pi\epsilon_0 r)$, becomes the associated Laguerre equation — another Sturm-Liouville problem. The completeness of the hydrogen eigenfunctions $R_{n\ell}(r)Y_\ell^m(\theta,\phi)$ is guaranteed by Sturm-Liouville theory, ensuring that any state of the hydrogen atom can be expanded in the energy eigenbasis.

5.4 Variational Characterization of Eigenvalues

The Rayleigh quotient for the Sturm-Liouville operator is:

$$\lambda[y] = \frac{\int_a^b \left[p(x)|y'|^2 + q(x)|y|^2\right]dx + \text{boundary terms}}{\int_a^b w(x)|y|^2\,dx}$$

The min-max theorem states that:

$$\lambda_1 = \min_{y \neq 0} \lambda[y], \qquad \lambda_n = \min_{\substack{y \neq 0 \\ y \perp y_1, \ldots, y_{n-1}}} \lambda[y]$$

In quantum mechanics, this is the variational principle: the ground state energy is the minimum of $\langle\psi|H|\psi\rangle / \langle\psi|\psi\rangle$ over all trial wave functions.

6. Singular Sturm-Liouville Problems

When $p(x)$ or $w(x)$ vanishes at an endpoint, or the interval is infinite, the problem is called singular. The boundary condition at a singular endpoint is replaced by requiring the solution to be bounded (or square-integrable).

Example: Legendre equation on $[-1, 1]$ with $p(x) = 1 - x^2$(vanishes at both endpoints):

$$-\frac{d}{dx}\!\left[(1-x^2)\frac{dy}{dx}\right] = \ell(\ell+1)\, y$$

Requiring boundedness at $x = \pm 1$ gives eigenvalues $\lambda_\ell = \ell(\ell+1)$with eigenfunctions $P_\ell(x)$. Despite the singular endpoints, all the SL theorems (real eigenvalues, orthogonality, completeness) still hold.

Example: Bessel equation on $[0, a]$ with $p(x) = x$ (vanishes at$x = 0$):

$$-\frac{d}{dx}\!\left(x\frac{dy}{dx}\right) + \frac{m^2}{x}\,y = \lambda\, x\, y$$

Requiring boundedness at $x = 0$ selects $J_m$ over $Y_m$. With$y(a) = 0$, the eigenvalues are $\lambda_{mn} = (j_{mn}/a)^2$.

7. Comparison Theorems and Asymptotics

7.1 Sturm Comparison Theorem

Theorem: Consider two equations $y'' + p_1(x)y = 0$ and$z'' + p_2(x)z = 0$ where $p_2(x) > p_1(x)$ on an interval. If$y$ has two consecutive zeros at $x_1$ and $x_2$, then $z$has at least one zero in $(x_1, x_2)$.

Physical interpretation: A stronger restoring force (larger $p_2$) leads to more rapid oscillation. In quantum mechanics, this tells us that deeper potential wells (lower $V$) produce more closely spaced nodes in the wave function.

7.2 Weyl's Asymptotic Law

For large eigenvalues, the density of eigenvalues follows a universal law. For the Dirichlet Laplacian on a domain $\Omega \subset \mathbb{R}^d$:

$$N(\lambda) \sim \frac{|\Omega|}{(4\pi)^{d/2}\,\Gamma(d/2 + 1)}\,\lambda^{d/2} \quad \text{as } \lambda \to \infty$$

where $N(\lambda)$ counts eigenvalues below $\lambda$ and$|\Omega|$ is the volume of the domain. This is Weyl's law, which famously answers the question "Can you hear the shape of a drum?" — the eigenvalue spectrum determines the area (and perimeter, via corrections) of the domain.

For the 1D problem $-y'' = \lambda y$ on $[0, L]$, Weyl's law gives$N(\lambda) \sim L\sqrt{\lambda}/\pi$, consistent with the exact eigenvalues$\lambda_n = (n\pi/L)^2$ giving $n \sim L\sqrt{\lambda}/\pi$.

7.3 WKB Approximation for Large Eigenvalues

For the Sturm-Liouville equation with large $\lambda$, the WKB (Wentzel-Kramers-Brillouin) approximation gives approximate eigenvalues via the quantization condition:

$$\int_a^b \sqrt{\lambda\, w(x)/p(x) - q(x)/p(x)}\,dx = \left(n - \frac{1}{2}\right)\pi + O(1/\sqrt{\lambda})$$

This is the Bohr-Sommerfeld quantization condition of old quantum theory, and provides remarkably accurate eigenvalue estimates even for moderate $n$.

7. Python Simulation: Sturm-Liouville Theory

This simulation numerically solves Sturm-Liouville eigenvalue problems, verifies orthogonality and completeness, demonstrates the Green's function expansion, and connects to quantum mechanics through the harmonic oscillator.

Python

script.py238 lines

#!/usr/bin/env python3
"""
sturm_liouville.py
Numerical exploration of Sturm-Liouville theory: eigenvalue problems,
orthogonality, completeness, Green's functions, and quantum mechanics.
numpy only (no scipy).
"""
import numpy as np

print("=" * 65)
print("  STURM-LIOUVILLE THEORY — NUMERICAL EXPLORATIONS")
print("=" * 65)
print()

# ─── Part 1: Numerical SL Eigenvalue Problem ────────────────────────
print("=== Part 1: Solving -y'' = lambda*y on [0, pi], y(0)=y(pi)=0 ===")
print()

N = 200  # number of grid points
a, b = 0, np.pi
h = (b - a) / (N + 1)
x = np.linspace(a + h, b - h, N)  # interior points

# Construct the finite-difference matrix for -d^2/dx^2
# -y'' approximated by (-y_{j-1} + 2y_j - y_{j+1}) / h^2
diag_main = 2.0 * np.ones(N) / h**2
diag_off = -1.0 * np.ones(N - 1) / h**2

# Build tridiagonal matrix
A = np.diag(diag_main) + np.diag(diag_off, 1) + np.diag(diag_off, -1)

# Solve eigenvalue problem
eigenvalues, eigenvectors = np.linalg.eigh(A)

print("Eigenvalues (first 10):")
print(f"  {'n':>3s}  {'Numerical':>14s}  {'Exact (n^2)':>14s}  {'Rel Error':>12s}")
print("  " + "-" * 48)
for n in range(1, 11):
    exact = n**2
    numerical = eigenvalues[n-1]
    rel_err = abs(numerical - exact) / exact
    print(f"  {n:3d}  {numerical:14.8f}  {exact:14.8f}  {rel_err:12.4e}")

print()

# ─── Part 2: Orthogonality Verification ─────────────────────────────
print("=== Part 2: Orthogonality of Eigenfunctions ===")
print()

# Normalize eigenvectors (with weight w=1)
norms = np.sqrt(np.sum(eigenvectors**2, axis=0) * h)
eigvecs_normalized = eigenvectors / norms[np.newaxis, :]

print("Inner products <y_n, y_m> (should be delta_nm):")
print(f"  {'(n,m)':>8s}  {'<y_n, y_m>':>14s}  {'Expected':>10s}")
print("  " + "-" * 36)
for n, m in [(1,1), (1,2), (2,2), (1,3), (3,3), (1,5), (5,5), (3,7)]:
    ip = np.sum(eigvecs_normalized[:, n-1] * eigvecs_normalized[:, m-1]) * h
    expected = 1.0 if n == m else 0.0
    print(f"  ({n:d},{m:d}){' '*(6-len(f'({n},{m})'))}{ip:14.8f}  {expected:10.1f}")

print()

# Compare with exact eigenfunctions
print("Comparison with exact eigenfunctions sin(n*pi*x/pi) = sin(n*x):")
for n in [1, 2, 3]:
    y_exact = np.sin(n * x)
    y_exact /= np.sqrt(np.sum(y_exact**2) * h)
    y_num = eigvecs_normalized[:, n-1]
    # Fix sign ambiguity
    if np.sum(y_num * y_exact) < 0:
        y_num = -y_num
    max_err = np.max(np.abs(y_num - y_exact))
    print(f"  n={n}: max|y_num - y_exact| = {max_err:.6e}")
print()

# ─── Part 3: Completeness (Eigenfunction Expansion) ─────────────────
print("=== Part 3: Completeness — Eigenfunction Expansion ===")
print()

# Expand f(x) = x*(pi-x) in eigenfunctions
f_target = x * (np.pi - x)

print("Expanding f(x) = x*(pi-x) in eigenfunctions:")
print(f"  {'N terms':>8s}  {'L2 error':>14s}  {'Max error':>14s}")
print("  " + "-" * 40)

for N_terms in [1, 3, 5, 10, 20, 50, 100]:
    N_use = min(N_terms, N)
    coeffs = np.zeros(N_use)
    for n in range(N_use):
        coeffs[n] = np.sum(eigvecs_normalized[:, n] * f_target) * h

f_approx = np.zeros(N)
    for n in range(N_use):
        f_approx += coeffs[n] * eigvecs_normalized[:, n]

l2_err = np.sqrt(np.sum((f_target - f_approx)**2) * h)
    max_err = np.max(np.abs(f_target - f_approx))
    print(f"  {N_terms:8d}  {l2_err:14.4e}  {max_err:14.4e}")

print()

# Show first few coefficients
print("First 6 expansion coefficients:")
for n in range(6):
    cn = np.sum(eigvecs_normalized[:, n] * f_target) * h
    # Exact: c_n = 4/(n*pi) for odd n (using sin series of x(pi-x))
    n_idx = n + 1
    exact_cn = 4.0 / (n_idx**3) if n_idx % 2 == 1 else 0.0
    # Note: normalization differs; adjust
    print(f"  c_{n_idx} = {cn:+.8f}")
print()

# ─── Part 4: Green's Function Expansion ─────────────────────────────
print("=== Part 4: Green's Function via Eigenfunction Expansion ===")
print()

# Solve -y'' = f(x) with y(0)=y(pi)=0 using Green's function
# f(x) = sin(x) => exact solution y = sin(x)/1 = sin(x) (lambda_1 = 1)
# More interesting: f(x) = 1 => y = x*(pi-x)/2

f_source = np.ones(N)  # f(x) = 1

# Eigenfunction expansion: y = sum c_n/lambda_n * y_n
y_green = np.zeros(N)
for n in range(min(100, N)):
    cn = np.sum(eigvecs_normalized[:, n] * f_source) * h
    y_green += cn / eigenvalues[n] * eigvecs_normalized[:, n]

# Exact solution: y = x*(pi-x)/2
y_exact_green = x * (np.pi - x) / 2

print("Solving -y'' = 1 via Green's function expansion:")
print(f"  Max error (100 modes): {np.max(np.abs(y_green - y_exact_green)):.4e}")
print(f"  y(pi/2) numerical:  {np.interp(np.pi/2, x, y_green):.8f}")
print(f"  y(pi/2) exact:      {(np.pi/2)*(np.pi - np.pi/2)/2:.8f}")
print()

# Convergence with number of modes
print("Convergence of Green's function expansion:")
for N_modes in [1, 5, 10, 25, 50, 100]:
    y_approx = np.zeros(N)
    for n in range(min(N_modes, N)):
        cn = np.sum(eigvecs_normalized[:, n] * f_source) * h
        y_approx += cn / eigenvalues[n] * eigvecs_normalized[:, n]
    err = np.max(np.abs(y_approx - y_exact_green))
    print(f"  {N_modes:3d} modes: max error = {err:.6e}")
print()

# ─── Part 5: Quantum Harmonic Oscillator ────────────────────────────
print("=== Part 5: Quantum Harmonic Oscillator (SL Problem) ===")
print()

# -psi'' + x^2*psi = E*psi (in units where hbar=m=omega=1, so E = 2n+1)
N_qm = 300
L_qm = 8.0  # large enough domain
h_qm = 2 * L_qm / (N_qm + 1)
x_qm = np.linspace(-L_qm + h_qm, L_qm - h_qm, N_qm)

# Kinetic energy: -d^2/dx^2
diag_main_qm = 2.0 / h_qm**2 + x_qm**2  # -d^2/dx^2 + V(x)
diag_off_qm = -1.0 * np.ones(N_qm - 1) / h_qm**2

H = np.diag(diag_main_qm) + np.diag(diag_off_qm, 1) + np.diag(diag_off_qm, -1)

E_vals, psi_vecs = np.linalg.eigh(H)

print("Energy eigenvalues (should be 2n+1 in natural units, i.e., 1, 3, 5, ...):")
print(f"  {'n':>3s}  {'E_numerical':>14s}  {'E_exact':>10s}  {'Error':>12s}")
print("  " + "-" * 44)
for n in range(8):
    E_exact = 2*n + 1
    E_num = E_vals[n]
    print(f"  {n:3d}  {E_num:14.8f}  {E_exact:10.4f}  {abs(E_num - E_exact):12.4e}")

print()

# Verify orthogonality
norms_qm = np.sqrt(np.sum(psi_vecs**2, axis=0) * h_qm)
psi_norm = psi_vecs / norms_qm[np.newaxis, :]

print("Orthogonality check <psi_n|psi_m>:")
for n, m in [(0,0), (0,1), (1,1), (0,2), (2,2), (1,3)]:
    ip = np.sum(psi_norm[:, n] * psi_norm[:, m]) * h_qm
    print(f"  <{n}|{m}> = {ip:+.8f}  (expected {1.0 if n==m else 0.0:+.1f})")
print()

# Node counting (oscillation theorem)
print("Node count (oscillation theorem: n-th state has n nodes):")
for n in range(6):
    psi_n = psi_norm[:, n]
    # Count zero crossings
    nodes = np.sum(np.diff(np.sign(psi_n)) != 0) // 2
    print(f"  n={n}: {nodes} nodes (expected {n})")
print()

# ─── Part 6: Rayleigh-Ritz Variational Principle ────────────────────
print("=== Part 6: Rayleigh Quotient (Variational Principle) ===")
print()

# For -y'' = lambda*y on [0, pi]
# Rayleigh quotient: R[y] = int y'^2 dx / int y^2 dx

print("Rayleigh quotient R[y] = int(y'^2)/int(y^2) for trial functions:")
print("  (minimum should be lambda_1 = 1.0)")
print()

# Trial 1: exact eigenfunction sin(x)
y_trial = np.sin(x)
dy = np.gradient(y_trial, h)
R = np.sum(dy**2) * h / (np.sum(y_trial**2) * h)
print(f"  sin(x):         R = {R:.8f}")

# Trial 2: parabola x*(pi-x) (not exact but close)
y_trial2 = x * (np.pi - x)
dy2 = np.gradient(y_trial2, h)
R2 = np.sum(dy2**2) * h / (np.sum(y_trial2**2) * h)
print(f"  x*(pi-x):       R = {R2:.8f}  (upper bound on lambda_1)")

# Trial 3: cubic x*(pi-x)^2
y_trial3 = x * (np.pi - x)**2
dy3 = np.gradient(y_trial3, h)
R3 = np.sum(dy3**2) * h / (np.sum(y_trial3**2) * h)
print(f"  x*(pi-x)^2:     R = {R3:.8f}")

# Trial 4: Gaussian (not satisfying BCs exactly, but close)
y_trial4 = np.exp(-(x - np.pi/2)**2 / 0.5)
y_trial4[0] = 0
y_trial4[-1] = 0
dy4 = np.gradient(y_trial4, h)
R4 = np.sum(dy4**2) * h / (np.sum(y_trial4**2) * h)
print(f"  Gaussian:       R = {R4:.8f}")

print()
print(f"  Exact lambda_1 = 1.0000000")
print(f"  All Rayleigh quotients >= lambda_1 (variational principle)")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

← Fourier Methods Part III Overview →