# What is the efficient way to code Successive Over-relaxation (SOR) method in Mathematica?

I have written a SOR method (SOR is this method) code using C-style procedural loops.

However, I think there might be (much) better ways to achieve the same end in Mma avoiding these loops.

So my question is: Is that feasible to code SOR method using functions like Nest, Fold that could be substantially more efficient?

For example, the following is a simple set of multiline equations with tridiagonal matrix.(my original equations will have an larger coefficient matrix, like 100*100 tridiagonal sparse array.) 1. Let listnew be the list that stores the unknown $x_1,x_2,x_3,x_4,x_5$ after each time of iteration.
2. And listold correspondingly stores the unknowns before each time of iteration.
3. b stores right-hand side of the equation above.
4. k indicates the times of iteration.
5. i is the index of $x_i$.
6. ϵ is the pre-defined maximum error.
7. error indicates the maximum difference between listnew and listold after each iteration.
8. ω is the so-called relaxation factor.

Then, my (dumb) code using For loop would be:

ϵ = 0.1;
ω = 1.5;
b = {1, 1, 1, 1, 1};
listnew = listold = Table[0, {i, 5}];
error = ϵ + 1;
For[k = 0, error > ϵ, k++,
listnew[] = listold[] - ω (b[] + 4 listold[] - listold[])/4;
For[i = 2, i < 5, i++,
listnew[[i]] = listold[[i]] -
ω (b[[i]] + 4 listold[[i]] - listnew[[i - 1]] - listold[[i + 1]])/4];
listnew[] = listold[] - ω (b[] + 4 listold[] - listnew[])/4;
error = Max[Abs[listnew - listold]]; listold = listnew;]


I'm sorry if it looks so dumb...

• Can you please explain what's SOR method and add your code? Apr 24 '14 at 5:45
• Can you share your implementation of SOR using the loop method in Mathematica? (may be with sample input, how it is called, and sample output?) SOR is: en.wikipedia.org/wiki/Successive_over-relaxation Apr 24 '14 at 6:04
• My response to this query has a simple implementation of the Jacobi method. It might be possible to modify that to get SOR. Apr 24 '14 at 14:08
• You are using for loop there over k, but k is not used anywhere in the body. You can replace For[k = 0, error > \[Epsilon], k++, by just While[error > \[Epsilon], Apr 25 '14 at 4:34
• If you just want to solve the equations rather than practicing programming SOR method, just use LinearSolve, your self-made solver will probably never have a chance to beat it. Apr 25 '14 at 4:53

You asking for "efficient" implementation. SOR is an iterative process. One does a step, checks if converged or not, and if not, repeat the step. I found M to be very fast using plain old Do. I do not see the point of replacing a well understand algorithm expressed naturally using Do[ f[i,....], {i,1,n}], by something like f[#,...] &/@ Range[n] or other variations, just for the sake of syntax change.

But the reason I am writing is not this. I just thought your iterative implementation was not clear so will post implementation of SOR based on direct implementation of the formula I am familiar with, as seen in this Mathematical association of America page The above solves A x = b using SOR. So $a_{ii}$ is the diagonal entry of $A$ and $k+1$ is the new solution and $k$ is the old solution. Guess starts with the zero solution.

The loop above is implemented using Do and those inner summations are implemented using the Sum command. The check for convergence uses relative error check. For additional speed, one can try to Compile this. As others said, SOR is not used in practice to solve real problem. It is only for academic interest as there are faster iterative methods these days, and if you have lots of RAM, you can just use a direct solver. LinearSolve in Mathematica.

Here is an implementation of the above equation. You can change omega and for fixed tolerance, see the effect on how many iterations are needed to converge. Manipulate[tick;
Module[{keepRunning = True, i, j, xOld = Table[0, {5}], n = 5},
While[keepRunning,
Do[x[[i]] = xOld[[i]] + w/mat[[i, i]] (b[[i]] -
Sum[mat[[i, j]] x[[j]], {j, 1, i - 1}] - Sum[mat[[i, j]] xOld[[j]], {j, i, n}])
, {i, 1, n}
];
If[Norm[x - xOld]/Norm[x] < e,
keepRunning = False,
xOld = x;
k++
]
]
];

Grid[{
{Row[{"Iteration needed to converge ", k}], SpanFromLeft},
{TableForm[{x, direct}, TableHeadings -> {{"SOR", "direct"},
{"x1", "x2", "x3", "x4", "x5"}}], SpanFromLeft}},
Spacings -> {1, 1}, Frame -> All],

Grid[{{Button["RUN", k = 0; x = Table[0, {5}]; tick = Not[tick]]}}],

{{w, 1.5, "w"}, 1.1, 1.9, .01, Appearance -> "Labeled"},(*sor omega*)
{{e, 0.0001, "error tolerance"}, 0.000000001, 0.001, 0.000000001,
Appearance -> "Labeled"},

{{x, Table[0, {5}]}, None},
{{k, 0}, None},
{{tick, False}, None},
TrackedSymbols :> {tick},

Initialization :> (SetOptions[\$FrontEndSession, PrintPrecision -> 16];
mat = {{-4, 1, 0, 0, 0}, {1, -4, 1, 0, 0}, {0, 1, -4, 1, 0}, {0, 0, 1, -4, 1},
{0, 0, 0, 1, -4}};
b = Table[1, {5}];
direct = N@LinearSolve[mat, b]
)
]

• Thanks for your thorough answer. Your answer clarified some of my confusions. And I found many to learn from your code. Based on my limited numerical analysis knowledge, I really don't know any modern (and more efficient) solver to implement in solving large matrix. I will for sure learn some more practical methods in the future. :-) Apr 27 '14 at 8:22

Nasser has given a good answer, here's just some complements.

First, a literal translation for the first formula provided in the wikipedia of SOR:

methodSOR1[a_, b_, ω_, ϕ_, ϵ_] :=
Module[{d, l, u, coe1, coe2, coe3}, {d, l, u} =
N@{DiagonalMatrix@Diagonal@#, LowerTriangularize[#, -1], UpperTriangularize[#, 1]} &@a;
coe1 = Inverse[d + ω l]; coe2 = ω b; coe3 = (ω u + (ω - 1) d);
FixedPoint[coe1.(coe2 - coe3.#) &, ϕ, SameTest -> (Max@Abs[#1 - #2] < ϵ &)]]


Then, the implementation of the second formula (compiled version), notice the Module inside the code is necessary for the speed up because it localizes intermediate variables i.e. ϕ, ϕold, n in the function:

methodSOR2 =
Compile[{{a, _Real, 2}, {b, _Real, 1}, {ϕ0, _Real, 1}, {ω, _Real}, {ϵ, _Real}},
Module[{ϕ = ϕ0, ϕold = ϕ0 + 1, n = Length@ϕ0},
While[Max@Abs[ϕold - ϕ] > ϵ, ϕold = ϕ;
Do[ϕ[[i]] = (1 - ω) ϕ[[i]] +
ω/a[[i, i]] (b[[i]] - Sum[a[[i, j]] ϕ[[j]], {j, i - 1}] -
Sum[a[[i, j]] ϕ[[j]], {j, i + 1, n}]), {i, n}]]; ϕ]];


Finally, since you have optimize the formula for your specific problem, I'd like to give an example of compiling your original code with minimal effort:

specificSOR =
With[{end = 5000},
With[{ϵ = 0.1, ω = 1.5, b = ConstantArray[1., {end}], list = ConstantArray[0., {end}]},
Compile[{},
Module[{listold = list, listnew = list, error = ϵ + 1, k, i},
For[k = 0, error > ϵ, k++,
listnew[] = listold[] - ω (b[] + 4 listold[] - listold[])/4;
For[i = 2, i < end, i++,
listnew[[i]] = listold[[i]] -
ω (b[[i]] + 4 listold[[i]] - listnew[[i - 1]] - listold[[i + 1]])/4];
listnew[[end]] = listold[[end]] -
ω (b[[end]] + 4 listold[[end]] - listnew[[end - 1]])/4;
error = Max[Abs[listnew - listold]]; listold = listnew];
listnew]]]];


Of course using Do and While instead of For will make the code look better.

Finally it bears repeating that if you're not practicing programming SOR, just use LinearSolve for your equations.