Solving the 1D Poisson Equation with Julia

This sample post solves a familiar boundary-value problem with a second-order finite-difference method. It also exercises the blog’s Markdown support: equations, Julia syntax highlighting, tables, blockquotes, lists, links, and inline code.

The model problem

Consider the one-dimensional Poisson equation

with homogeneous Dirichlet boundary conditions $u(0)=u(1)=0$. The exact solution is

On a uniform grid with spacing $h$, the centered approximation is

Because the truncation error is $O(h^2)$, halving the grid spacing should reduce the error by approximately a factor of four.

Julia implementation

The matrix is symmetric and tridiagonal, so Julia’s SymTridiagonal type stores and solves the system efficiently.

using LinearAlgebra

function solve_poisson_1d(interior_points)
    grid_spacing = 1 / (interior_points + 1)
    x = collect(range(
        grid_spacing,
        1 - grid_spacing;
        length=interior_points,
    ))

    operator = SymTridiagonal(
        fill(2 / grid_spacing^2, interior_points),
        fill(-1 / grid_spacing^2, interior_points - 1),
    )
    right_hand_side = pi^2 .* sin.(pi .* x)
    numerical_solution = operator \ right_hand_side
    exact_solution = sin.(pi .* x)

    return maximum(abs.(numerical_solution .- exact_solution))
end

for interior_points in (32, 64, 128, 256)
    error = solve_poisson_1d(interior_points)
    println("n = $interior_points, max error = $error")
end

The code uses only Julia’s standard library. A run with Julia produced:

n = 32,  max error = 7.547363359678982e-4
n = 64,  max error = 1.9463264608221564e-4
n = 128, max error = 4.9421936692550794e-5
n = 256, max error = 1.2452237266202815e-5

Convergence results

The observed order approaches two, which is the expected behavior. This gives us three quick checks:

1. The boundary conditions are represented by excluding the two boundary points from the unknowns.
2. The numerical solution converges to the known analytical solution.
3. The convergence rate agrees with the stencil’s $O(h^2)$ truncation error.

This is a deliberately small example, but the same workflow scales to more interesting operators: define the discretization, exploit its structure, and verify it against a result you know. Browse the Julia tag for related posts.¹