Solve a 2D Poisson equation problem using Gauss-Seidel

gauss_seidel.go

// Solve the Poisson equation for a 2D grid using the Gauss-Seidel method
// Gauss-Seidel method: https://en.wikipedia.org/wiki/Gauss–Seidel_method

// The problem used is Example 1 from https://my.ece.utah.edu/~ece6340/LECTURES/Feb1/Nagel%202012%20-%20Solving%20the%20Generalized%20Poisson%20Equation%20using%20FDM.pdf
// A 4x4 grid of electrical values
//   - top boundary is Dirichlet boundary fixed at 1.0V
//   - bottom boundary is Dirichlet boundary fixed at 0.0V
//   - left & right boundaries are Neumann boundaries fixed at derivative 0.0V

package main

import (
	"fmt"
	"math"
)

// Transform vector x by matrix A
func transform(A [][]float64, x []float64) []float64 {
	newX := make([]float64, len(A))
	for i := range A {
		for j := range A[i] {
			newX[i] += A[i][j] * x[j]
		}
	}
	return newX
}

// Calculate the residual: a scalar measure for how close the
// transformed guess x is to the right hand side b. The larger
// the residual the further apart.
func residual(A [][]float64, x, b []float64) float64 {
	var residual float64

	newX := transform(A, x)
	for i := range newX {
		delta := (b[i] - newX[i])
		residual += delta * delta
	}

	return residual
}

func main() {
	// Build the coefficient matrix A
	// [ 1  0  0  0  0  0  0  0  0  0  0  0 ]
	// [ 0  1  0  0  0  0  0  0  0  0  0  0 ]
	// [ 0  0  1 -1  0  0  0  0  0  0  0  0 ]
	// [ 1  0  1 -4  1  0  0  1  0  0  0  0 ]
	// [ 0  1  0  1 -4  1  0  0  1  0  0  0 ]
	// [ 0  0  0  0 -1  1  0  0  0  0  0  0 ]
	// [ 0  0  0  0  0  0  1 -1  0  0  0  0 ]
	// [ 0  0  0  1  0  0  1 -4  1  0  1  0 ]
	// [ 0  0  0  0  1  0  0  1 -4  1  0  1 ]
	// [ 0  0  0  0  0  0  0  0 -1  1  0  0 ]
	// [ 0  0  0  0  0  0  0  0  0  0  1  0 ]
	// [ 0  0  0  0  0  0  0  0  0  0  0  1 ]
	A := make([][]float64, 12)
	for i := range A {
		A[i] = make([]float64, 12)
	}
	A[0][0] = 1
	A[1][1] = 1
	A[2][2] = 1
	A[2][3] = -1
	A[3][0] = 1
	A[3][2] = 1
	A[3][3] = -4
	A[3][4] = 1
	A[3][7] = 1
	A[4][1] = 1
	A[4][3] = 1
	A[4][4] = -4
	A[4][5] = 1
	A[4][8] = 1
	A[5][4] = -1
	A[5][5] = 1
	A[6][6] = 1
	A[6][7] = -1
	A[7][3] = 1
	A[7][6] = 1
	A[7][7] = -4
	A[7][8] = 1
	A[7][10] = 1
	A[8][4] = 1
	A[8][7] = 1
	A[8][8] = -4
	A[8][9] = 1
	A[8][11] = 1
	A[9][8] = -1
	A[9][9] = 1
	A[10][10] = 1
	A[11][11] = 1

	// Build the b vector
	// [0 0 0 0 0 0 0 0 0 0 1 1]^T
	b := make([]float64, 12)
	b[10] = 1
	b[11] = 1

	// Use Gauss-Seidel to compute the unknown vector x
	// A . x = b
	// x = A^-1 . b
	x := make([]float64, 12)
	res := math.MaxFloat64
	iter := 0
	for ; iter < 1000 && res > 1e-5; iter++ {
		for i := 0; i < len(A); i++ {
			var sum float64

			for j := 0; j < len(A[i]); j++ {
				if j != i {
					sum += A[i][j] * x[j]
				}
			}
			x[i] = (b[i] - sum) / A[i][i]
		}

		res = residual(A, x, b)
	}
	fmt.Printf("Iter: %d, Residual: %f\n", iter, res)
	fmt.Printf("Final: %v\n", x)
	fmt.Printf("Actual: %v\n", []float64{0, 0, 1.0 / 3, 1.0 / 3, 1.0 / 3, 1.0 / 3, 2.0 / 3, 2.0 / 3, 2.0 / 3, 2.0 / 3, 1, 1})
}

Iter: 15, Residual: 0.000007
Final: [0 0 0.3304064869880676 0.3314176835119724 0.3317835386842489 0.3317835386842489 0.6642987877130508 0.6651168707758188 0.66541285533458 0.66541285533458 1 1]
Actual: [0 0 0.3333333333333333 0.3333333333333333 0.3333333333333333 0.3333333333333333 0.6666666666666666 0.6666666666666666 0.6666666666666666 0.6666666666666666 1 1]