PDE Cookbook
A collection of worked examples solving classical partial differential equations with findiff. Each recipe is self-contained.
1D Forced Harmonic Oscillator
Solve
\[u'' - \alpha\, u' + \omega^2\, u = \cos(2t)\]
on \([0, 10]\) with \(u(0)=0,\; u(10)=1\).
import numpy as np
from findiff import Diff, Identity, PDE, BoundaryConditions
N = 300
t = np.linspace(0, 10, N)
dt = t[1] - t[0]
alpha = 1.0
omega = np.sqrt(5)
L = Diff(0, dt)**2 - alpha * Diff(0, dt) + omega**2 * Identity()
f = np.cos(2 * t)
bc = BoundaryConditions((N,))
bc[0] = 0
bc[-1] = 1
pde = PDE(L, f, bc)
u = pde.solve()
2D Laplace Equation (Steady-State Heat)
Solve \(\nabla^2 u = 0\) on a unit square with temperature fixed on all four edges:
import numpy as np
from findiff import Diff, PDE, BoundaryConditions
N = 80
x = y = np.linspace(0, 1, N)
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, :] = 0 # left edge: u = 0
bc[-1, :] = 0 # right edge: u = 0
bc[:, 0] = 0 # bottom edge: u = 0
bc[:, -1] = np.sin(np.pi * x) # top edge: u = sin(pi*x)
pde = PDE(L, f, bc)
u = pde.solve()
2D Poisson Equation
Solve \(\nabla^2 u = -2\pi^2 \sin(\pi x)\sin(\pi y)\) on a unit square with \(u=0\) on all boundaries (exact solution: \(u = \sin(\pi x)\sin(\pi y)\)).
import numpy as np
from findiff import Diff, PDE, BoundaryConditions
N = 80
x = y = np.linspace(0, 1, N)
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 = -2 * np.pi**2 * np.sin(np.pi * X) * np.sin(np.pi * Y)
bc = BoundaryConditions(X.shape)
bc[0, :] = 0
bc[-1, :] = 0
bc[:, 0] = 0
bc[:, -1] = 0
pde = PDE(L, f, bc)
u = pde.solve()
# Check
u_exact = np.sin(np.pi * X) * np.sin(np.pi * Y)
print("Max error:", np.max(np.abs(u - u_exact)))
2D Heat Conduction with Mixed BCs
A plate with a temperature profile on one edge, zero heat flux across the others:
import numpy as np
from findiff import Diff, PDE, BoundaryConditions
N = 100
x = np.linspace(0, 1, N)
y = np.linspace(0, 1, N)
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[1, :] = Diff(0, dx), 0 # Neumann: du/dx = 0 at x = 0
bc[-1, :] = 300. - 200 * Y # Dirichlet at x = 1
bc[:, 0] = 300. # Dirichlet at y = 0
bc[1:-1, -1] = Diff(1, dy), 0 # Neumann: du/dy = 0 at y = 1
pde = PDE(L, f, bc)
u = pde.solve()
1D Schrodinger Eigenvalue Problem
Find the lowest energy levels of the quantum harmonic oscillator \(-\tfrac{1}{2}\psi'' + \tfrac{1}{2}x^2\psi = E\psi\):
import numpy as np
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)
print("Eigenvalues:", eigenvalues)
# Expected: ≈ [0.5, 1.5, 2.5, 3.5, 4.5, 5.5]
Advection-Diffusion with Robin BC
Solve
\[D\, u'' - v\, u' = 0\]
on \([0, 1]\) with \(u(0)=1\) and a Robin condition \(u + \tfrac{D}{v} u' = 0\) at \(x=1\) (outflow condition):
import numpy as np
from findiff import Diff, PDE, BoundaryConditions
N = 200
x = np.linspace(0, 1, N)
dx = x[1] - x[0]
D = 0.01 # diffusion coefficient
v = 1.0 # advection velocity
L = D * Diff(0, dx)**2 - v * Diff(0, dx)
f = np.zeros(N)
bc = BoundaryConditions((N,))
bc[0] = 1 # Dirichlet: u(0) = 1
bc[-1] = (1, Diff(0, dx), D / v, 0) # Robin: u + (D/v)*u' = 0
pde = PDE(L, f, bc)
u = pde.solve()
Large 2D Problem with Iterative Solver
For large grids, iterative solvers with preconditioning can be faster and use less memory than the default direct solver:
import numpy as np
from findiff import Diff, PDE, BoundaryConditions
N = 300
x = y = np.linspace(0, 1, N)
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.ones_like(X)
bc = BoundaryConditions(X.shape)
bc[0, :] = 0
bc[-1, :] = 0
bc[:, 0] = 0
bc[:, -1] = 0
pde = PDE(L, f, bc)
u = pde.solve(solver='cg', preconditioner='ilu', tol=1e-10)
1D Heat Equation (Time-Dependent)
Solve \(u_t = D\, u_{xx}\) on \([0, 1]\) with \(u_0 = \sin(\pi x)\) and homogeneous Dirichlet conditions:
import numpy as np
from findiff import Diff, TimeDependentPDE, BoundaryConditions
nx = 101
x = np.linspace(0, 1, nx)
dx = x[1] - x[0]
D = 0.01
L = D * Diff(0, dx)**2
u0 = np.sin(np.pi * x)
bc = BoundaryConditions((nx,))
bc[0] = 0
bc[-1] = 0
t = np.linspace(0, 1, 500)
u_final = TimeDependentPDE(L, u0, bc, t).solve()
u_exact = np.sin(np.pi * x) * np.exp(-D * np.pi**2 * t[-1])
print("Max error:", np.max(np.abs(u_final - u_exact)))
2D Diffusion (Time-Dependent)
Solve \(u_t = D\,(\nabla^2 u)\) on a unit square with zero boundary conditions:
import numpy as np
from findiff import Diff, TimeDependentPDE, BoundaryConditions
nx, ny = 51, 51
x = np.linspace(0, 1, nx)
y = np.linspace(0, 1, ny)
dx, dy = x[1] - x[0], y[1] - y[0]
X, Y = np.meshgrid(x, y, indexing='ij')
D = 0.01
L = D * (Diff(0, dx)**2 + Diff(1, dy)**2)
u0 = np.sin(np.pi * X) * np.sin(np.pi * Y)
bc = BoundaryConditions((nx, ny))
bc[0, :] = 0
bc[-1, :] = 0
bc[:, 0] = 0
bc[:, -1] = 0
t = np.linspace(0, 0.5, 200)
u_final = TimeDependentPDE(L, u0, bc, t).solve()
u_exact = (np.sin(np.pi * X) * np.sin(np.pi * Y)
* np.exp(-2 * D * np.pi**2 * t[-1]))
print("Max error:", np.max(np.abs(u_final - u_exact)))