The Min–Max theorem¶

As described in the 18.303 class notes, the min–max theorem (or "variational theorem") says that, for a self-adjoint operator $\hat{A}$ under some inner product $\langle u, v \rangle$ on some Sobolev space $u,v \in V$ (essentially, $u,v$ such that the "bilinear form" $\langle u, \hat{A} v \rangle$ is finite), the Rayleigh quotient $$R\{u\} = \frac{\langle u, \hat{A} u \rangle}{\langle u, u \rangle}$$ evaluated for any $u \ne 0$ in $V$ is an upper bound for the minimum eigenvalue of $\hat{A}$ (if any) and a lower bound for the maximum eigenvalue of $\hat{A}$ (if any). Furthermore, the eigenvalues of $\hat{A}$ are extrema of $R\{u\}$.

As described in the notes, this allows us to easily guess the qualitative features of the first few eigenvalues of an operator, even in cases that we can't easily solve analytically.

Finite-difference discretization¶

As an example, let's consider $\hat{A} = -\nabla\cdot c\nabla$ for some $c(\vec{x}) > 0$, which is self-adjoint and positive-definite for any finite domain $\Omega$ with Dirichlet boundary conditions $u|_{\partial\Omega}=0$ and the usual inner product $\langle u, v \rangle = \int_\Omega \bar{u} v$.

We'll define the finite-difference discretization by code that is a bit more general than we have used before: we'll allow the user to pass in 1d arrays x and y of the $x$ and $y$ coordinates used in the discretization (possibly nonuniformly spaced!), including the coordinates of the endpoints where $u$ is zero. We'll also allow the user to pass in an arbitrary function $c(x,y)$ that we'll evaluate on the grid as needed.

In [2]:
# compute the first-derivative finite-difference matrix
# for Dirichlet boundaries, given a grid x[:] of x points
# (including the endpoints where the function = 0!).
function sdiff1(x)
N = length(x) - 2
dx1 = Float64[1/(x[i+1] - x[i]) for i = 1:N]
dx2 = Float64[-1/(x[i+1] - x[i]) for i = 2:N+1]
spdiagm((dx1,dx2), (0,-1), N+1, N)
end

flatten(X) = reshape(X, length(X))

# compute the -∇⋅ c ∇ operator for a function c(x,y)
# and arrays x[:] and y[:] of the x and y points,
# including the endpoints where functions are zero
# (i.e. Dirichlet boundary conditions).
function Laplacian(x, y, c)
Dx = sdiff1(x)
Nx = size(Dx,2)
Dy = sdiff1(y)
Ny = size(Dy,2)

G = [kron(speye(Ny), Dx); kron(Dy, speye(Nx))]

# grids for derivatives in x and y directions
x′ = [0.5*(x[i]+x[i+1]) for i = 1:length(x)-1]
y′ = [0.5*(y[i]+y[i+1]) for i = 1:length(y)-1]

# evaluate c(x)
C = spdiagm([ flatten(Float64[c(X,Y) for X in x′, Y in y[2:end-1]]);
flatten(Float64[c(X,Y) for X in x[2:end-1], Y in y′]) ])

return G' * C * G # -∇⋅ c ∇
end

Out[2]:
Laplacian (generic function with 1 method)

As a quick check to make sure this is working, let's do $\hat{A}=-\nabla^2$ on the box $[0,1]\times[0,1]$, and verify that it gives the expected separable solution and first eigenvalue $2\pi^2$:

In [9]:
x = linspace(0,1,100)
A = Laplacian(x, x, (x,y) -> 1.0);
λ, X = eigs(A, nev=1, which=:SM);

using PyPlot
U = reshape(X[:,1],98,98)
pcolor(x[2:end]',x[2:end],U)
axis("equal")