The Castleford projectRelaxation in Python: Jacobi, Gauss-Seidel and over-relaxation
Hide part of a real height map, refill it by averaging, and measure how fast and how well each method works.
The grid comes from the OpenTopoData service, which serves the Copernicus EU-DEM elevation model. We asked for 61 by 61 points about 50 m apart, centred on Castleford, and got back 3,721 heights between 5.6 m and 74.7 m. The learner then cuts nine square holes, each 19 points by 19, so that 361 heights at a time disappear, and tries to put them back using nothing but the ring of known heights around each hole.
The rule for every missing point is the same: become the average of your four neighbours. There is exactly one set of numbers that satisfies it everywhere at once, and Python can find it directly with a sparse matrix solver, which gives us a reference answer to measure against. The interesting part is getting there by sweeping. Jacobi's method updates every point from the old values of its neighbours. Gauss-Seidel uses each new value the moment it is computed. Successive over-relaxation, or SOR, does the same but deliberately overshoots each update by a factor, here 1.729, the figure theory recommends for a hole of this size.
Sweeps needed across the nine hidden blocks (our Python, EU-DEM heights)| Method | Sweeps to get within 1 cm of the exact fill | Where "stop when nothing moves 1 cm" stops | How far off it still is then |
|---|
| Jacobi | 219 to 356 | 45 to 135 | 0.27 to 0.43 m |
| Gauss-Seidel | 111 to 214 | 33 to 84 | 0.16 to 0.30 m |
| SOR, factor 1.729 | 24 to 40 | 23 to 41 | no more than 0.01 m |
Two lessons come out of that table. The first is speed. On average Gauss-Seidel needed about half the sweeps of Jacobi, and SOR needed about a fifth of Gauss-Seidel's. One line of code, multiplying each correction by 1.729, cut the work by that much. The second lesson is subtler and more useful. The natural stopping rule, quit when no point moves by more than a centimetre, is badly wrong for the slow methods. Jacobi's updates shrink to a centimetre while its answer is still up to 43 cm from where it will end up; each step is small only because the method is slow, not because it has arrived. SOR is the only one of the three for which small steps really do mean done.
How close the refilled heights come to the real ones (our comparison)| Fill | Error per block, root mean square | Closer of the two |
|---|
| Relaxation from the ring | 1.25 to 4.74 m | 9 blocks of 9 |
| Every hole set to the ring's average | 1.96 to 16.71 m | 0 blocks of 9 |
Relaxation always beat the lazy fill, but it has a built-in blind spot. An average can never be larger than the largest thing it averages, so the refilled patch can never rise above the highest point of its ring or sink below the lowest. In two of the nine blocks the real ground inside rose above every point on the ring, and relaxation smoothed that high ground away. The worst single height was off by 16.3 m. No amount of sweeping fixes that; the information was simply not in the ring.
Ages 8 to 11
Fill a paper grid with missing numbers by averaging neighbours, round after round, and watch the numbers settle.
Ages 11 to 15
Write the Jacobi sweep with two nested loops in Python and count the rounds on a small grid.
Ages 15 and up
Use NumPy on the real heights, add Gauss-Seidel and SOR, and design a stopping rule that tells the truth.
Data and caveats
Heights from the EU-DEM v1.1 model through the OpenTopoData API, produced using Copernicus data and information funded by the European Union. A 25 m elevation model smooths buildings and banks, and our points are about 50 m apart, so none of these figures describes a particular street. The holes, the fills and every count above come from our own code.
