1D Poisson by Finite Differences

Discretize a boundary-value problem into a tridiagonal system and solve for the potential across the interval.

Problem

The 1D Poisson equation -u''(x) = f(x) with fixed boundary values is the simplest elliptic boundary-value problem. Finite differences turn the second derivative into a tridiagonal linear system that we solve directly, the building block for many field solvers.

Discretization

On a uniform grid with spacing h, -u''(x_i) is approximated by (-u_{i-1} + 2 u_i - u_{i+1})/h^2. This produces a matrix with 2 on the diagonal and -1 on the off-diagonals, scaled by 1/h^2. Boundary values move to the right-hand side.

import numpy as np
n=6; h=1/(n+1); x=np.linspace(h,1-h,n)
f=np.ones(n)                      # source -u'' = 1
A=(2*np.eye(n)-np.eye(n,k=1)-np.eye(n,k=-1))/h**2
u=np.linalg.solve(A,f)
exact=0.5*x*(1-x)                 # analytic solution
print('max error',round(np.max(np.abs(u-exact)),6))

Result

With f=1 and zero boundaries the exact solution is the parabola u(x)=x(1-x)/2, and the finite-difference answer matches it to machine precision because the scheme is exact for quadratics. For general sources the error shrinks as h^2 (second order). The tridiagonal system is solved in linear time by the Thomas algorithm, so even fine grids are cheap in 1D.

• The discrete Laplacian is symmetric positive definite, making iterative solvers like CG applicable in higher dimensions.
• Second-order accuracy means halving h quarters the error.
• This operator, in 2D and 3D form, is the core of the electrostatic and magnetostatic solves in Kronos plasma models.