pde/laplace_2d.viv

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