Solving Partial Differential Equations
findiff can solve linear partial differential equations with Dirichlet, Neumann or Robin boundary conditions using sparse matrix representations.
Setup
import numpy as np
from findiff import Diff, PDE, BoundaryConditions
The PDE class takes a differential operator and solves the resulting
linear system on a given grid.
Example: 1D Boundary Value Problem
Solve \(u''(x) = \sin(x)\) on \([0, \pi]\) with \(u(0) = 0\) and \(u(\pi) = 0\):
x = np.linspace(0, np.pi, 100)
dx = x[1] - x[0]
# The differential equation: u'' = sin(x)
L = Diff(0, dx)**2
f = np.sin(x)
# Boundary conditions
bc = BoundaryConditions(x.shape)
bc[0] = 0 # u(0) = 0 (Dirichlet)
bc[-1] = 0 # u(pi) = 0 (Dirichlet)
# Solve
pde = PDE(L, f, bc)
u = pde.solve()
Example: 2D Laplace Equation
Solve \(\nabla^2 u = 0\) on a square domain:
x = y = np.linspace(0, 1, 50)
dx, dy = x[1] - x[0], y[1] - y[0]
X, Y = np.meshgrid(x, y, indexing='ij')
L = Diff(0, dx)**2 + Diff(1, dy)**2
f = np.zeros_like(X)
bc = BoundaryConditions(X.shape)
bc[0, :] = X[0, :] # Dirichlet at x = 0
bc[-1, :] = X[-1, :] # Dirichlet at x = 1
bc[:, 0] = Y[:, 0] # Dirichlet at y = 0
bc[:, -1] = Y[:, -1] # Dirichlet at y = 1
pde = PDE(L, f, bc)
u = pde.solve()
Neumann Boundary Conditions
For Neumann BCs, pass a Diff operator as the second element:
bc = BoundaryConditions(x.shape)
bc[0] = 0 # u(0) = 0
bc[-1] = Diff(0, dx), 1.0 # u'(L) = 1.0
Robin (Mixed) Boundary Conditions
Robin boundary conditions have the form
\(\alpha \, u + \beta \, \frac{\partial u}{\partial n} = g\).
Specify them as a 4-tuple (alpha, diff_op, beta, g):
bc = BoundaryConditions(x.shape)
bc[0] = 1 # Dirichlet: u(0) = 1
bc[-1] = (1, Diff(0, dx), 1, 3) # Robin: u + u' = 3
For 2D problems, Robin BCs on an edge work the same way:
bc = BoundaryConditions(X.shape)
bc[-1, :] = (1, Diff(0, dx), 0.5, g_boundary) # alpha*u + beta*du/dx = g
bc[0, :] = u_left # Dirichlet on other edges
bc[:, 0] = u_bottom
bc[:, -1] = u_top
Alternatively, you can construct the Robin operator yourself using
Identity() and Diff and pass it as a 2-tuple (like a Neumann BC):
robin_op = alpha * Identity() + beta * Diff(0, dx)
bc[-1] = robin_op, g
Eigenvalue Problems
Differential operators can solve eigenvalue problems of the form
\(L\,u = \lambda\,u\) using the eigsh method (for symmetric
operators) or eigs (for general operators). These are thin wrappers
around scipy.sparse.linalg.eigsh / eigs.
Example: vibration modes of a string
The eigenvalue problem \(u'' = \lambda\,u\) on \([0, \pi]\) with \(u(0) = u(\pi) = 0\) has exact eigenvalues \(\lambda_n = -n^2\):
x = np.linspace(0, np.pi, 200)
dx = x[1] - x[0]
L = Diff(0, dx)**2
bc = BoundaryConditions(x.shape)
bc[0] = 0
bc[-1] = 0
eigenvalues, eigenvectors = L.eigsh(x.shape, k=5, which='SM', bc=bc)
# eigenvalues ≈ [-1, -4, -9, -16, -25]
Pass bc to eliminate boundary degrees of freedom (homogeneous
Dirichlet). Eigenvectors are returned with shape (*grid_shape, k)
— access the i-th mode as eigenvectors[..., i].
For generalized eigenvalue problems \(L\,u = \lambda\,M\,u\), pass a
second operator via the M parameter:
eigenvalues, eigenvectors = L.eigsh(shape, k=5, which='SM', bc=bc, M=M_op)
See the Matrix Representation guide for a more detailed example with the Schrodinger equation.
Iterative Solvers
By default, PDE.solve() uses the direct solver
scipy.sparse.linalg.spsolve. For large problems (especially in 3D),
iterative solvers can be more memory-efficient and faster. Pass the
solver keyword to select an iterative method:
u = pde.solve(solver='gmres')
The following solver names are supported:
Name |
Method |
When to use |
|---|---|---|
|
|
Default. Reliable for small to medium problems. |
|
Conjugate Gradient |
Symmetric positive definite systems only. |
|
Generalized Minimal Residual |
General (non-symmetric) systems. Good default iterative choice. |
|
BiConjugate Gradient Stabilized |
General non-symmetric systems. |
|
LGMRES |
Variant of GMRES with improved convergence in some cases. |
|
Minimal Residual |
Symmetric indefinite systems. |
Solver options such as tolerance, maximum iterations, and initial guess are passed as additional keyword arguments:
u = pde.solve(solver='gmres', tol=1e-10, maxiter=1000)
# Provide an initial guess close to the expected solution
u = pde.solve(solver='bicgstab', x0=u_previous)
Preconditioners can dramatically improve convergence. Use the
preconditioner='ilu' shorthand for an incomplete LU factorization:
u = pde.solve(solver='gmres', preconditioner='ilu')
Or pass a custom scipy.sparse.linalg.LinearOperator:
from scipy.sparse.linalg import LinearOperator
M = LinearOperator(...) # your custom preconditioner
u = pde.solve(solver='gmres', preconditioner=M)
Custom solver callables are also supported. Any function with
signature f(A, b, **kw) -> x or f(A, b, **kw) -> (x, info)
can be used:
from scipy.sparse.linalg import spsolve
u = pde.solve(solver=lambda A, b: spsolve(A, b))
If an iterative solver fails to converge, a RuntimeError is raised.