Skip to content

Commit d80c031

Browse files
committed
Prettified SOR source
1 parent 2ca5603 commit d80c031

2 files changed

Lines changed: 32 additions & 28 deletions

File tree

Successive_Over_Relaxation/Lua/Yonaba/sor.lua

Lines changed: 12 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -4,26 +4,27 @@
44
-- Creates a vector of values
55
local function vector(len, v)
66
local x = {}
7-
for i = 1,len do x[i] = v end
7+
for i = 1, len do x[i] = v end
88
return x
99
end
1010

1111
-- Computes the norm of a given vector
1212
local function norm(v)
1313
local n = 0
14-
for i,_v in ipairs(v) do
14+
for i, _v in ipairs(v) do
1515
n = n + (_v * _v)
1616
end
1717
return math.sqrt(n)
1818
end
1919

2020
-- Solves the given matrix using SOR algorithm
21-
-- m : a matrix
21+
-- m : a matrix A in [A]*[x] = [b]
22+
-- b : the result vector in [A]*[x] = [b]
2223
-- w : the relaxation parameter, defaults to 1.86
2324
-- eps : the convergence criterion, defaults to 1e-6
2425
-- maxiter : the maximum number of iterations, defaults to 1e4
2526
-- return : a solution vector
26-
local function SOR(m, w, eps, maxiter)
27+
local function SOR(m, b, w, eps, maxiter)
2728
local n = #m
2829
local x = vector(n, 0)
2930
local q, p, sum
@@ -36,17 +37,17 @@ local function SOR(m, w, eps, maxiter)
3637
repeat
3738
t = t + 1
3839
q = norm(x)
39-
for i=1,n do
40-
sum = m[i][n+1]
41-
for j=1,n do
42-
if (i~=j) then
43-
sum = sum - m[i][j]*x[j]
40+
for i =1, n do
41+
sum = 0
42+
for j = 1, n do
43+
if (i ~= j) then
44+
sum = sum + m[i][j] * x[j]
4445
end
4546
end
46-
x[i] = (1-w) * x[i] + w * sum/m[i][i]
47+
x[i] = (1 - w) * x[i] + (w / m[i][i]) * (b[i] - sum)
4748
end
4849
p = norm(x)
49-
until (math.abs(p-q) < eps) or (t >= maxiter)
50+
until (math.abs(p - q) < eps) or (t >= maxiter)
5051

5152
return x
5253
end

Successive_Over_Relaxation/Lua/Yonaba/sor_test.lua

Lines changed: 20 additions & 17 deletions
Original file line numberDiff line numberDiff line change
@@ -17,32 +17,35 @@ local function run(message, f)
1717
print(('%02d. %68s: %s'):format(total, dec(message,68), status))
1818
end
1919

20+
-- Fuzzy equality test
2021
local function fuzzy_eq(a, b, eps)
2122
eps = eps or 1e-6
2223
return math.abs(a - b) < eps
2324
end
2425

2526
run('Testing SOR', function()
2627
local m = {
27-
{-4,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,-11},
28-
{1,-4,1,0,0,1,0,0,0,0,0,0,0,0,0,0,-03},
29-
{0,1,-4,1,0,0,1,0,0,0,0,0,0,0,0,0,-03},
30-
{0,0,1,-4,0,0,0,1,0,0,0,0,0,0,0,0,-11},
31-
{1,0,0,0,-4,1,0,0,1,0,0,0,0,0,0,0,-08},
32-
{0,1,0,0,1,-4,1,0,0,1,0,0,0,0,0,0, 00},
33-
{0,0,1,0,0,1,-4,1,0,0,1,0,0,0,0,0, 00},
34-
{0,0,0,1,0,0,1,-4,0,0,0,1,0,0,0,0,-08},
35-
{0,0,0,0,1,0,0,0,-4,1,0,0,1,0,0,0,-08},
36-
{0,0,0,0,0,1,0,0,1,-4,1,0,0,1,0,0, 00},
37-
{0,0,0,0,0,0,1,0,0,1,-4,1,0,0,1,0, 00},
38-
{0,0,0,0,0,0,0,1,0,0,1,-4,0,0,0,1,-08},
39-
{0,0,0,0,0,0,0,0,1,0,0,0,-4,1,0,0,-10},
40-
{0,0,0,0,0,0,0,0,0,1,0,0,1,-4,1,0,-02},
41-
{0,0,0,0,0,0,0,0,0,0,1,0,0,1,-4,1,-02},
42-
{0,0,0,0,0,0,0,0,0,0,0,1,0,0,1,-4,-10},
28+
{-4,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0},
29+
{1,-4,1,0,0,1,0,0,0,0,0,0,0,0,0,0},
30+
{0,1,-4,1,0,0,1,0,0,0,0,0,0,0,0,0},
31+
{0,0,1,-4,0,0,0,1,0,0,0,0,0,0,0,0},
32+
{1,0,0,0,-4,1,0,0,1,0,0,0,0,0,0,0},
33+
{0,1,0,0,1,-4,1,0,0,1,0,0,0,0,0,0},
34+
{0,0,1,0,0,1,-4,1,0,0,1,0,0,0,0,0},
35+
{0,0,0,1,0,0,1,-4,0,0,0,1,0,0,0,0},
36+
{0,0,0,0,1,0,0,0,-4,1,0,0,1,0,0,0},
37+
{0,0,0,0,0,1,0,0,1,-4,1,0,0,1,0,0},
38+
{0,0,0,0,0,0,1,0,0,1,-4,1,0,0,1,0},
39+
{0,0,0,0,0,0,0,1,0,0,1,-4,0,0,0,1},
40+
{0,0,0,0,0,0,0,0,1,0,0,0,-4,1,0,0},
41+
{0,0,0,0,0,0,0,0,0,1,0,0,1,-4,1,0},
42+
{0,0,0,0,0,0,0,0,0,0,1,0,0,1,-4,1},
43+
{0,0,0,0,0,0,0,0,0,0,0,1,0,0,1,-4},
4344
}
4445

45-
local x = sor(m)
46+
local b = {-11, -3, -3, -11, -8, 0, 0, -8, -8, 0, 0, -8, -10,-2, -2, -10}
47+
48+
local x = sor(m, b)
4649
assert(fuzzy_eq(x[1], 5.454459))
4750
assert(fuzzy_eq(x[2], 4.594688))
4851
assert(fuzzy_eq(x[3], 4.594679))

0 commit comments

Comments
 (0)