Added Gauss Seidel in Lua

Yonaba · Yonaba · commit 79ba7914e4df · 2014-03-25T17:36:45.000Z
diff --git a/Gauss_Seidel_Algorithm/Lua/Yonaba/gauss_seidel.lua b/Gauss_Seidel_Algorithm/Lua/Yonaba/gauss_seidel.lua
@@ -0,0 +1,55 @@
+-- Gauss Seidel algorithm implementation
+-- See : http://en.wikipedia.org/wiki/Gauss%E2%80%93Seidel_method
+
+-- Creates a vector of values
+local function vector(len, v)
+  local x = {}
+  for i = 1, len do x[i] = v end
+  return x
+end
+
+-- Computes the norm of a given vector
+local function norm(v)
+  local n = 0
+  for i, _v in ipairs(v) do
+    n = n + (_v * _v)
+  end
+  return math.sqrt(n)
+end
+
+-- Solves the given matrix using Gauss-Seidel algorithm
+-- m       : a matrix A in [A][x] = [b]
+-- b       : the result vector in [A][x] = [b]
+-- w       : the relaxation parameter, defaults to 1.86
+-- eps     : the convergence criterion, defaults to 1e-6
+-- maxiter : the maximum number of iterations, defaults to 1e4
+-- return  : a solution vector
+local function gauss_seidel(m, b, w, eps, maxiter)
+  local n = #m
+  local x = vector(n, 0)
+  local q, p, sum
+  local t = 0
+
+  w = w or 1.86
+  eps = eps or 1e-6
+  maxiter = maxiter or 1e4
+
+  repeat
+    t = t + 1
+    q = norm(x)
+    for i = 1, n do
+      sum = 0
+      for j = 1, n do
+        if (i ~= j) then
+          sum = sum + m[i][j] * x[j]
+        end
+      end
+      x[i] = (1 / m[i][i]) * (b[i] - sum)
+    end
+    p = norm(x)
+  until (math.abs(p - q) < eps) or (t >= maxiter)
+
+  return x
+end
+
+return gauss_seidel
diff --git a/Gauss_Seidel_Algorithm/Lua/Yonaba/gauss_seidel_test.lua b/Gauss_Seidel_Algorithm/Lua/Yonaba/gauss_seidel_test.lua
@@ -0,0 +1,69 @@
+-- Tests for sor.lua
+local gauss_seidel = require 'gauss_seidel'
+
+local total, pass = 0, 0
+
+local function dec(str, len)
+  return #str < len
+     and str .. (('.'):rep(len-#str))
+      or str:sub(1,len)
+end
+
+local function run(message, f)
+  total = total + 1
+  local ok, err = pcall(f)
+  if ok then pass = pass + 1 end
+  local status = ok and 'PASSED' or 'FAILED'
+  print(('%02d. %68s: %s'):format(total, dec(message,68), status))
+end
+
+-- Fuzzy equality test
+local function fuzzy_eq(a, b, eps)
+  eps = eps or 1e-3
+  return math.abs(a - b) < eps
+end
+
+run('Testing Gauss Seidel', function()
+  local m = {
+    {-4,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0},
+    {1,-4,1,0,0,1,0,0,0,0,0,0,0,0,0,0},
+    {0,1,-4,1,0,0,1,0,0,0,0,0,0,0,0,0},
+    {0,0,1,-4,0,0,0,1,0,0,0,0,0,0,0,0},
+    {1,0,0,0,-4,1,0,0,1,0,0,0,0,0,0,0},
+    {0,1,0,0,1,-4,1,0,0,1,0,0,0,0,0,0},
+    {0,0,1,0,0,1,-4,1,0,0,1,0,0,0,0,0},
+    {0,0,0,1,0,0,1,-4,0,0,0,1,0,0,0,0},
+    {0,0,0,0,1,0,0,0,-4,1,0,0,1,0,0,0},
+    {0,0,0,0,0,1,0,0,1,-4,1,0,0,1,0,0},
+    {0,0,0,0,0,0,1,0,0,1,-4,1,0,0,1,0},
+    {0,0,0,0,0,0,0,1,0,0,1,-4,0,0,0,1},
+    {0,0,0,0,0,0,0,0,1,0,0,0,-4,1,0,0},
+    {0,0,0,0,0,0,0,0,0,1,0,0,1,-4,1,0},
+    {0,0,0,0,0,0,0,0,0,0,1,0,0,1,-4,1},
+    {0,0,0,0,0,0,0,0,0,0,0,1,0,0,1,-4},
+  }
+
+  local b = {-11, -3, -3, -11, -8, 0, 0, -8, -8, 0, 0, -8, -10,-2, -2, -10}
+
+  local x = gauss_seidel(m, b)
+  assert(fuzzy_eq(x[1], 5.454459))
+  assert(fuzzy_eq(x[2], 4.594688))
+  assert(fuzzy_eq(x[3], 4.594679))
+  assert(fuzzy_eq(x[4], 5.454531))
+  assert(fuzzy_eq(x[5], 6.223469))
+  assert(fuzzy_eq(x[6], 5.329580))
+  assert(fuzzy_eq(x[7], 5.329561))
+  assert(fuzzy_eq(x[8], 6.223491))
+  assert(fuzzy_eq(x[9], 6.109833))
+  assert(fuzzy_eq(x[10], 5.170488))
+  assert(fuzzy_eq(x[11], 5.170455))
+  assert(fuzzy_eq(x[12], 6.109845))
+  assert(fuzzy_eq(x[13], 5.045460))
+  assert(fuzzy_eq(x[14], 4.071980))
+  assert(fuzzy_eq(x[15], 4.071964))
+  assert(fuzzy_eq(x[16], 5.045457))
+end)
+
+print(('-'):rep(80))
+print(('Total : %02d: Pass: %02d - Failed : %02d - Success: %.2f %%')
+  :format(total, pass, total-pass, (pass*100/total)))
diff --git a/Gauss_Seidel_Algorithm/tags b/Gauss_Seidel_Algorithm/tags
@@ -0,0 +1 @@
+Algebra, Numerical, Matrix, Linear Systems

Original file line number	Diff line number	Diff line change
`@@ -0,0 +1 @@`
	`1`	`+Algebra, Numerical, Matrix, Linear Systems`