PDE solvers · Steady states and modes (elliptic)

Laplace's equation

\nabla^{2} u = 0

Method: successive over-relaxation, with the optimal relaxation factor for a square
Checked: a value symmetry fixes exactly; and that over-relaxation beats plain Gauss-Seidel by the factor theory predicts

A square plate: one edge held at 1, the other three at 0. What is the temperature at the center?

Exactly $\frac{1}{4}$ , and no series is needed to see it. Solve the problem four times, heating each edge in turn, and add the answers: the sum has every edge at 1, so it is 1 everywhere. By symmetry each of the four gives the same value at the center. The discrete equations have the same symmetry, so the numerical answer is $\frac{1}{4}$ too — to the solver's tolerance, not the grid's. It is the rare test with no discretization error in it.

Gauss-Seidel replaces each point by the average of its neighbors, sweep after sweep, and is slow: the error shrinks by a factor of about $1 - π^{2} h^{2}$ per sweep. SOR moves each point a factor $ω$ past that average. With $ω = 2/ (1 + sin π h)$ the factor becomes about $1 - 2 π h$ , and the number of sweeps drops from $O (n^{2})$ to $O (n)$ .

What the run shows

A plate with one hot edge

25 × 25 grid

the center of the plate is exactly 1/4, by symmetry
Gauss-Seidel reaches the same answer
optimal over-relaxation needs an order of magnitude fewer sweeps

OMEGA: 1.76909
center: 0.25
sor_sweeps: 110
gs_sweeps: 1315
speedup: 11.9545

Computed when this page was built, by the same interpreter that runs it in your browser: 3 of 3 tests, 2.13 s. The numbers above are the program's show lines.

The program

63 lines, one function. It is editable: change it and press Run (or ⌘/Ctrl + Enter). The tests at the bottom run too, so if your change breaks the solver, a test fails.

relax: Relax to convergence with relaxation factor omega; returns the grid.

Viv [1]

pi = 3.141592653589793
SWEEPS = zeros(1)       #: SWEEPS  how many sweeps the last solve took

#| Relax to convergence with relaxation factor omega; returns the grid.
def relax(n, omega, tol) {
    u = zeros(n + 2.0, n + 2.0)
    for j to n + 2.0 {
        u[0.0, j] = 1.0
    }
    change = 1.0
    sweeps = 0.0
    while change > tol and sweeps < 20000.0 {
        change = 0.0
        for i to n {
            for j to n {
                a = i + 1.0
                b = j + 1.0
                avg = 0.25 * (u[a - 1.0, b] + u[a + 1.0, b] + u[a, b - 1.0] + u[a, b + 1.0])
                step = omega * (avg - u[a, b])
                u[a, b] = u[a, b] + step
                change = max(change, abs(step))
            }
        }
        sweeps = sweeps + 1.0
    }
    SWEEPS[0] = sweeps
    return u
}

N = 23
H = 1.0 / (N + 1.0)
OMEGA = 2.0 / (1.0 + sin(pi * H))
u_sor = relax(N, OMEGA, 1.0e-12)
sor_sweeps = SWEEPS[0]
u_gs = relax(N, 1.0, 1.0e-12)
gs_sweeps = SWEEPS[0]
center = u_sor[(N + 1.0) / 2.0, (N + 1.0) / 2.0]
speedup = gs_sweeps / sor_sweeps
show OMEGA
show center
show sor_sweeps
show gs_sweeps
show speedup

# transposed, the first index of u_sor runs across and the second up, so the hot
# edge, u[0, j] = 1, is on the left
figure("A plate with one hot edge", "x", "y")
origin("lower")
plot2d(transpose(u_sor), "temperature")

test "the center of the plate is exactly 1/4, by symmetry" {
    assert_near(center, 0.25, 1.0e-10)
}
test "Gauss-Seidel reaches the same answer" {
    assert_near(u_gs[(N + 1.0) / 2.0, (N + 1.0) / 2.0], 0.25, 1.0e-9)
}
test "optimal over-relaxation needs an order of magnitude fewer sweeps" {
    assert_gt(speedup, 10.0)
    assert_lt(sor_sweeps, 5.0 * N)
}

# SPDX-License-Identifier: MIT OR Apache-2.0
# SPDX-FileCopyrightText: 2026 The Viv authors

pi = 3.141592653589793
SWEEPS = zeros(1)       #: SWEEPS  how many sweeps the last solve took

#| Relax to convergence with relaxation factor omega; returns the grid.
def relax(n, omega, tol) {
    u = zeros(n + 2.0, n + 2.0)
    for j to n + 2.0 {
        u[0.0, j] = 1.0
    }
    change = 1.0
    sweeps = 0.0
    while change > tol and sweeps < 20000.0 {
        change = 0.0
        for i to n {
            for j to n {
                a = i + 1.0
                b = j + 1.0
                avg = 0.25 * (u[a - 1.0, b] + u[a + 1.0, b] + u[a, b - 1.0] + u[a, b + 1.0])
                step = omega * (avg - u[a, b])
                u[a, b] = u[a, b] + step
                change = max(change, abs(step))
            }
        }
        sweeps = sweeps + 1.0
    }
    SWEEPS[0] = sweeps
    return u
}

N = 23
H = 1.0 / (N + 1.0)
OMEGA = 2.0 / (1.0 + sin(pi * H))
u_sor = relax(N, OMEGA, 1.0e-12)
sor_sweeps = SWEEPS[0]
u_gs = relax(N, 1.0, 1.0e-12)
gs_sweeps = SWEEPS[0]
center = u_sor[(N + 1.0) / 2.0, (N + 1.0) / 2.0]
speedup = gs_sweeps / sor_sweeps
show OMEGA
show center
show sor_sweeps
show gs_sweeps
show speedup

# transposed, the first index of u_sor runs across and the second up, so the hot
# edge, u[0, j] = 1, is on the left
figure("A plate with one hot edge", "x", "y")
origin("lower")
plot2d(transpose(u_sor), "temperature")

test "the center of the plate is exactly 1/4, by symmetry" {
    assert_near(center, 0.25, 1.0e-10)
}
test "Gauss-Seidel reaches the same answer" {
    assert_near(u_gs[(N + 1.0) / 2.0, (N + 1.0) / 2.0], 0.25, 1.0e-9)
}
test "optimal over-relaxation needs an order of magnitude fewer sweeps" {
    assert_gt(speedup, 10.0)
    assert_lt(sor_sweeps, 5.0 * N)
}

# SPDX-License-Identifier: MIT OR Apache-2.0
# SPDX-FileCopyrightText: 2026 The Viv authors

Open in the workbench The atlas Download laplace_2d.viv