Matrix Representation

Diff objects can be converted to sparse matrix representations. This is useful for eigenvalue problems or when you need to solve linear systems involving differential operators.

Basic usage

import numpy as np
from findiff import Diff

x = np.linspace(0, 6, 7)
d2_dx2 = Diff(0, x[1] - x[0]) ** 2

mat = d2_dx2.matrix(x.shape)
print(mat.toarray())

Output:

[[ 2. -5.  4. -1.  0.  0.  0.]
 [ 1. -2.  1.  0.  0.  0.  0.]
 [ 0.  1. -2.  1.  0.  0.  0.]
 [ 0.  0.  1. -2.  1.  0.  0.]
 [ 0.  0.  0.  1. -2.  1.  0.]
 [ 0.  0.  0.  0.  1. -2.  1.]
 [ 0.  0.  0. -1.  4. -5.  2.]]

The returned matrix is a scipy.sparse matrix.

Eigenvalue Problems

Differential operators have built-in eigs and eigsh methods for solving eigenvalue problems directly.

Example: 1D quantum harmonic oscillator

\[-\frac{1}{2}\frac{d^2\psi}{dx^2} + \frac{1}{2}x^2\psi = E\,\psi\]

The exact eigenvalues are \(E_n = n + \tfrac{1}{2}\).

from findiff import Diff, Identity, BoundaryConditions

n = 300
x = np.linspace(-8, 8, n)
dx = x[1] - x[0]

H = -0.5 * Diff(0, dx)**2 + 0.5 * x**2 * Identity()

bc = BoundaryConditions((n,))
bc[0] = 0
bc[-1] = 0

eigenvalues, eigenvectors = H.eigsh((n,), k=6, which='SM', bc=bc)
# eigenvalues ≈ [0.5, 1.5, 2.5, 3.5, 4.5, 5.5]

The eigsh method is for symmetric operators (Laplacian, Hamiltonian). Use eigs for general (non-symmetric) operators.

Key parameters:

k: number of eigenvalues to compute
which: which eigenvalues — 'SM' (smallest magnitude), 'SA' (smallest algebraic), 'LM' (largest magnitude)
bc: a BoundaryConditions object; boundary DOFs are eliminated (homogeneous Dirichlet)
sigma: shift for shift-invert mode (finds eigenvalues near sigma)
M: another operator for generalized problems \(L\psi = \lambda M\psi\)

Eigenvectors are returned with shape (*grid_shape, k), so eigenvectors[..., i] is the i-th eigenvector on the grid.

Example: 3D Schrodinger equation

\[-\frac{1}{2}\nabla^2 \psi + V\,\psi = E\,\psi\]

laplace = Diff(0, dx)**2 + Diff(1, dy)**2 + Diff(2, dz)**2
H = -0.5 * laplace + V * Identity()

bc = BoundaryConditions(shape)
bc[0, :, :] = 0
bc[-1, :, :] = 0
bc[:, 0, :] = 0
bc[:, -1, :] = 0
bc[:, :, 0] = 0
bc[:, :, -1] = 0

eigenvalues, eigenvectors = H.eigsh(shape, k=6, which='SM', bc=bc)