SOR Method

Click on the link for the interactive version:

Consider the matrix equation $$A\mathbf{x}=\mathbf{b}.$$ where $$A=\left( \begin{array}{cccc} 1&0&0&0\\ -1&2&-1&0\\ 0&-1&2&-1\\ 0&0&-1&1 \end{array} \right) \text{, } \mathbf{x}=\left( \begin{array}{c} x_0\\ x_1\\ x_2\\ x3 \end{array} \right)\text{ and } \mathbf{b}=\left( \begin{array}{c} 1\\ \frac14\\ \frac12\\ 0 \end{array} \right) .$$ Open the example code on my website click here, it will provide a template to implement the code outlined in this demonstration. The final solution is also on the web click here. Out aim is to implement SOR to find the solution $\mathbf{x}$. The task will be split into several steps:

1. Copy the code into a new project, check it compiles and runs.
2. Update your guess for $\mathbf{x}$ by rearranging each line of the matrix equation to get an equation for $x_0$, $x_1$, $x_2$ and $x_3$. For example $$x_0=(b_0-A_{0,1}x_1-A_{0,2}x_2-A_{0,3}x_3)/A_{0,0}.$$ Put the four expressions for each value of $x_j$ into the code at the appropriate place, take advantage of the zeros in the matrix to simplify formula.
3. Run the code, does the solution converge?
4. Try outputing the error by summing up the residuals. The residual is given by $$r_i=b_i-\sum_jA_{i,j}x_j.$$ Can you see the residual going down as your guess approaches the solution?
5. Now try to use the parameter ωω to overrelax your new guess with the formula $$x^{q+1}_j=x^q_j+\omega(x^{q+1}_j-x^q_j).$$ Can you choose $\omega$ to speed up convergence?

To solve this problem first we need to declare the appropriate libraries at the top of the code:

#include <iostream>
#include <cmath>
#include <vector>
using namespace std;

Inside the main function we need to set up the problem, we use three vectors a, b, c for the diagonals of the matrix, and a vector rhs for $\mathbf{b}$.

// declare vector to store matrix A
vector<double> a = { 0.,-1,-1,-1};
vector<double> b = { 1.,2,2,1};
vector<double> c = { 0.,-1,-1,0.};
// declare vector for the rhs v and solution x 
vector<double> rhs = {1.,0.25,0.5,0.};
vector<double> x;

Now we set up the matrix solver using a simple iteration scheme. The general code layout is outlined below, at the moment the code breaks after one step.

{
    // initial guess for x
    x = {0.,0.,0.,0.};
    // MATRIX SOLVE
    // sor variables
    // set iterMax to a large enough value so that method converges
    int iterMax=50;
    double tol=1.e-6,omega=1;
    // sor loop
    for(int sor=0;sor<iterMax;sor++)
    {
        // implement sor in here
        // x[0]=...
    
        // output guess at (sor+1)th iteration
        cout << "x = { ";
        for(auto xi : x)cout << xi << " ";
        cout << "} at iteration " << sor+1 << "\n";
        // make an exit condition when solution found
        // ...
        break;
    }
}

x = { 0 0 0 0 } at iteration 1

Now input the sor code, exploiting the tridiagonal structure of the matrix so simplify the process. We can see that the solution converges after about 50 steps (to 6sf).

{
    // initial guess for x
    x = {0.,0.,0.,0.};
    // MATRIX SOLVE
    // sor variables
    // set iterMax to a large enough value so that method converges
    int iterMax=50;
    double tol=1.e-6,omega=1;
    int sor=0;
    // sor loop
    for(sor=0;sor<iterMax;sor++)
    {
        // implement sor in here
        x[0] = rhs[0]/b[0];
        for(int j=1;j<3;j++)
        {
            x[j] = (rhs[j] - a[j]*x[j-1] - c[j]*x[j+1])/b[j];
        }
        x[3] = (rhs[3] - a[3]*x[2])/b[3];
        // output guess at (sor+1)th iteration
        cout << "x = { ";
        for(auto xi : x)cout << xi << " ";
        cout << "} at iteration " << sor+1 << "\n";
    }
}

x = { 1 0.625 0.5625 0.5625 } at iteration 1
x = { 1 0.90625 0.984375 0.984375 } at iteration 2
x = { 1 1.11719 1.30078 1.30078 } at iteration 3
x = { 1 1.27539 1.53809 1.53809 } at iteration 4
x = { 1 1.39404 1.71606 1.71606 } at iteration 5
x = { 1 1.48303 1.84955 1.84955 } at iteration 6
x = { 1 1.54977 1.94966 1.94966 } at iteration 7
x = { 1 1.59983 2.02475 2.02475 } at iteration 8
x = { 1 1.63737 2.08106 2.08106 } at iteration 9
x = { 1 1.66553 2.12329 2.12329 } at iteration 10
x = { 1 1.68665 2.15497 2.15497 } at iteration 11
x = { 1 1.70249 2.17873 2.17873 } at iteration 12
x = { 1 1.71436 2.19655 2.19655 } at iteration 13
x = { 1 1.72327 2.20991 2.20991 } at iteration 14
x = { 1 1.72995 2.21993 2.21993 } at iteration 15
x = { 1 1.73497 2.22745 2.22745 } at iteration 16
x = { 1 1.73872 2.23309 2.23309 } at iteration 17
x = { 1 1.74154 2.23732 2.23732 } at iteration 18
x = { 1 1.74366 2.24049 2.24049 } at iteration 19
x = { 1 1.74524 2.24286 2.24286 } at iteration 20
x = { 1 1.74643 2.24465 2.24465 } at iteration 21
x = { 1 1.74732 2.24599 2.24599 } at iteration 22
x = { 1 1.74799 2.24699 2.24699 } at iteration 23
x = { 1 1.74849 2.24774 2.24774 } at iteration 24
x = { 1 1.74887 2.24831 2.24831 } at iteration 25
x = { 1 1.74915 2.24873 2.24873 } at iteration 26
x = { 1 1.74937 2.24905 2.24905 } at iteration 27
x = { 1 1.74952 2.24929 2.24929 } at iteration 28
x = { 1 1.74964 2.24946 2.24946 } at iteration 29
x = { 1 1.74973 2.2496 2.2496 } at iteration 30
x = { 1 1.7498 2.2497 2.2497 } at iteration 31
x = { 1 1.74985 2.24977 2.24977 } at iteration 32
x = { 1 1.74989 2.24983 2.24983 } at iteration 33
x = { 1 1.74992 2.24987 2.24987 } at iteration 34
x = { 1 1.74994 2.2499 2.2499 } at iteration 35
x = { 1 1.74995 2.24993 2.24993 } at iteration 36
x = { 1 1.74996 2.24995 2.24995 } at iteration 37
x = { 1 1.74997 2.24996 2.24996 } at iteration 38
x = { 1 1.74998 2.24997 2.24997 } at iteration 39
x = { 1 1.74998 2.24998 2.24998 } at iteration 40
x = { 1 1.74999 2.24998 2.24998 } at iteration 41
x = { 1 1.74999 2.24999 2.24999 } at iteration 42
x = { 1 1.74999 2.24999 2.24999 } at iteration 43
x = { 1 1.75 2.24999 2.24999 } at iteration 44
x = { 1 1.75 2.24999 2.24999 } at iteration 45
x = { 1 1.75 2.25 2.25 } at iteration 46
x = { 1 1.75 2.25 2.25 } at iteration 47
x = { 1 1.75 2.25 2.25 } at iteration 48
x = { 1 1.75 2.25 2.25 } at iteration 49
x = { 1 1.75 2.25 2.25 } at iteration 50

Next we want to know how to track convergence of our method. Often one of the best ways is to looks at a residual analysis, which is basically where we calculate $$\mathbf{r} = \mathbf{b} - A\mathbf{x} .$$

We have to do this after we have solved for all elements in the vector $\mathbf{x}$.

{
    // initial guess for x
    x = {0.,0.,0.,0.};
    // MATRIX SOLVE
    // sor variables
    // set iterMax to a large enough value so that method converges
    int iterMax=50;
    double tol=1.e-6,omega=1;
    int sor=0;
    // sor loop
    for(sor=0;sor<iterMax;sor++)
    {
        vector<double> r(4);
        // implement sor in here
        x[0] = rhs[0]/b[0];
        for(int j=1;j<3;j++)
        {
            x[j] = (rhs[j] - a[j]*x[j-1] - c[j]*x[j+1])/b[j];
        }
        x[3] = (rhs[3] - a[3]*x[2])/b[3];
        // output guess at (sor+1)th iteration
        if(sor%10==9)
        {
            cout << "x = { ";
            for(auto xi : x)cout << xi << " ";
            cout << "} and r = {";
            cout << rhs[0] - b[0]*x[0] << " ";
            cout << rhs[1] - a[1]*x[0] - b[1]*x[1] - c[1]*x[2] << " ";
            cout << rhs[2] - a[2]*x[1] - b[2]*x[2] - c[2]*x[3] << " ";
            cout << rhs[3] - a[3]*x[2] - b[3]*x[3] << " ";
            cout << "} at iteration " << sor+1 << "\n";
        }
    }
}

x = { 1 1.66553 2.12329 2.12329 } and r = {0 0.0422351 0.0422351 0 } at iteration 10
x = { 1 1.74524 2.24286 2.24286 } and r = {0 0.00237841 0.00237841 0 } at iteration 20
x = { 1 1.74973 2.2496 2.2496 } and r = {0 0.000133937 0.000133937 0 } at iteration 30
x = { 1 1.74998 2.24998 2.24998 } and r = {0 7.54244e-06 7.54244e-06 0 } at iteration 40
x = { 1 1.75 2.25 2.25 } and r = {0 4.24741e-07 4.24741e-07 0 } at iteration 50

Here we can see that the size of the residual is broadly similar to the error in the solution. As a result we can use the size of the residual (i.e. the norm) to measure how close we are to the solution. We introduce an error term to the calculation and check it against the tolerance. Further we also implement SOR, storing our new estimate in a variable y so that we can generate a new value x[j] = x[j] + omega*(y-x[j]).

{
    // initial guess for x
    x = {0.,0.,0.,0.};
  // MATRIX SOLVE
  // sor variables
  // set iterMax to a large enough value so that method converges
  int iterMax=500;
  double tol=1.e-6,omega=1.4;
  int sor=0;
  // sor loop
  for(sor=0;sor<iterMax;sor++)
  {
    double error=0.;
    // implement sor in here
    {
      double y = rhs[0]/b[0];
      x[0] = x[0] + omega*(y-x[0]); 
    }
    for(int j=1;j<3;j++)
    {
      double y = (rhs[j] - a[j]*x[j-1] - c[j]*x[j+1])/b[j];
      x[j] = x[j] + omega*(y-x[j]); 
    }
    {
      double y = (rhs[3] - a[3]*x[2])/b[3];
      x[3] = x[3] + omega*(y-x[3]); 
    }
    // calculate residual norm ||r|| as sum of absolute values
    error += fabs(rhs[0] - b[0]*x[0]);
    error += fabs(rhs[1] - a[1]*x[0] - b[1]*x[1] - c[1]*x[2]);
    error += fabs(rhs[2] - a[2]*x[1] - b[2]*x[2] - c[2]*x[3]);
    error += fabs(rhs[3] - a[3]*x[2] - b[3]*x[3]);
    // make an exit condition when solution found
    if(error<tol)
      break;
  }
  cout << "Converged after " << sor+1 << " iterations.\n";
}
x

Converged after 19 iterations.

{ 1, 1.75, 2.25, 2.25 }

We can easily convert this into a function, which we can put in a separate file click here or above the main function. The function is called sorSolve(...) where the arguments are a, b, c, rhs, x, iterMax, tol, omega and sor. The definition of each parameter can be implied from the code above.

/* 
ON INPUT:
a, b and c -- are the diagonals of the matrix
rhs        -- is the right hand side
x          -- is the initial guess
iterMax    -- is maximum iterations
tol        -- is the tolerance level
omega      -- is the relaxation parameter 
sor        -- not used
ON OUTPUT:
a, b, c, rhs        -- unchanged
x                   -- solution to Ax=b
iterMax, tol, omega -- unchanged
sor                 -- number of iterations to converge
*/
void sorSolve(const std::vector<double> &a,const std::vector<double> &b,const std::vector<double> &c,const std::vector<double> &rhs,
std::vector<double> &x,int iterMax,double tol,double omega,int &sor )
{
    // assumes vectors a,b,c,d,rhs and x are same size (doesn't check)
    int n=a.size()-1;
    // sor loop
  for(sor=0;sor<iterMax;sor++)
  {
    double error=0.;
    // implement sor in here
    {
      double y = (rhs[0] - c[0]*x[1])/b[0];
      x[0] = x[0] + omega*(y-x[0]); 
    }
    for(int j=1;j<n;j++)
    {
      double y = (rhs[j] - a[j]*x[j-1] - c[j]*x[j+1])/b[j];
      x[j] = x[j] + omega*(y-x[j]); 
    }
    {
      double y = (rhs[n] - a[n]*x[n-1])/b[n];
      x[n] = x[n] + omega*(y-x[n]); 
    }
    // calculate residual norm ||r|| as sum of absolute values
    error += std::fabs(rhs[0] - b[0]*x[0] - c[0]*x[1]);
    for(int j=1;j<n;j++) 
        error += std::fabs(rhs[j] - a[j]*x[j-1] - b[j]*x[j] - c[j]*x[j+1]);
    error += std::fabs(rhs[n] - a[n]*x[n-1] - b[n]*x[n]);
    // make an exit condition when solution found
    if(error<tol)
      break;
  }
}

Once in a function we can play around to see the effect of the different input parameters. For instance, with defaults as above we get:

// reset initial guess for x
x = {0.,0.,0.,0.};
int sor;
// iterMax=50, tol=1.e-6, omega = 1
sorSolve(a,b,c,rhs,x,50,1.e-6,1.,sor);
cout << "Converged after " << sor+1 << " iterations.\n";
x

Converged after 50 iterations.

{ 1, 1.75, 2.25, 2.25 }

So with omega=1 we converge in 50 iterations. Now try again with omega=1.1:

// reset initial guess for x
x = {0.,0.,0.,0.};
// this time omega = 1.1
sorSolve(a,b,c,rhs,x,50,1.e-6,1.1,sor);
cout << "Converged after " << sor+1 << " iterations.\n";
x

Converged after 40 iterations.

{ 1, 1.75, 2.25, 2.25 }

Now try it with a few different values

for(int varyOmega=5;varyOmega<=20;varyOmega++)
{
    // reset initial guess for x
    x = {0.,0.,0.,0.};
    // this time omega = 1.1
    double omega = varyOmega/10.;
    sorSolve(a,b,c,rhs,x,5000,1.e-6,omega,sor);
    cout << "Try with omega = " << omega << " we get x = { ";
    for(auto xi : x)cout << xi << " ";
    cout << "} after " << sor+1 << " iterations.\n";
}

Try with omega = 0.5 we get x = { 1 1.75 2.25 2.25 } after 158 iterations.
Try with omega = 0.6 we get x = { 1 1.75 2.25 2.25 } after 122 iterations.
Try with omega = 0.7 we get x = { 1 1.75 2.25 2.25 } after 97 iterations.
Try with omega = 0.8 we get x = { 1 1.75 2.25 2.25 } after 78 iterations.
Try with omega = 0.9 we get x = { 1 1.75 2.25 2.25 } after 62 iterations.
Try with omega = 1 we get x = { 1 1.75 2.25 2.25 } after 50 iterations.
Try with omega = 1.1 we get x = { 1 1.75 2.25 2.25 } after 40 iterations.
Try with omega = 1.2 we get x = { 1 1.75 2.25 2.25 } after 31 iterations.
Try with omega = 1.3 we get x = { 1 1.75 2.25 2.25 } after 21 iterations.
Try with omega = 1.4 we get x = { 1 1.75 2.25 2.25 } after 19 iterations.
Try with omega = 1.5 we get x = { 1 1.75 2.25 2.25 } after 26 iterations.
Try with omega = 1.6 we get x = { 1 1.75 2.25 2.25 } after 36 iterations.
Try with omega = 1.7 we get x = { 1 1.75 2.25 2.25 } after 53 iterations.
Try with omega = 1.8 we get x = { 1 1.75 2.25 2.25 } after 87 iterations.
Try with omega = 1.9 we get x = { 1 1.75 2.25 2.25 } after 192 iterations.
Try with omega = 2 we get x = { 0 -3331.25 4 -3329.5 } after 5001 iterations.

As we can see the number of iterations is very dependent on the choice of omega. If we choose the value $\omega=1.4$ the solver will converge around 2.5 times quicker, whilst choosing a value that is too big of $\omega=1.9$ takes around 4 times longer. In fact, theoretically there is always an optimal choice of $\omega\in[1,2)$ that depend on the properties of the matrix $A$ (the spectral radius). For some simple matrices this can be calculated easily but in general it is difficult to know.