Successive Over Relaxation (SOR) method

It looks like we had a longer session in the
last lecture, so will try to make it slightly lesser this time. To begin with we are looking
at advance method which will increase the rate of convergence of our iteration method. When we say increase over what we can get
with Jacobi method and Gauss-Seidel method and one of the first methods that were proposed
is what is known as the successive over relaxation method. This method is, goes together with Jacobi
method and Gauss-Seidel method as almost like a basic iterative method, because it is not
very different from those methods. So, the idea of this successive over relaxation is
that, in the general case in an iterative method we have some way of determining x k
plus 1 from x k and there is an increment, we are in essence incrementing in x k by delta
x k to get x k plus 1. We know that in the case of convergent schemes and especially
under asymptotic convergence conditions and we are approaching this solution asymptotically
from one side. So, if you are approaching asymptotically
then that means, that we are taking steps, but may be not taking as big as step as can
be taken and we are taking small steps. So, if you can increase the step size more than
what we think it should be, because we are getting delta x k from the equation and the
amount of delta x k that we get depends on the iteration matrix p and it also depend
on x k. So, that means that delta x k is not constant, the step size is not constant, but
we seem to be every time under estimating that step size. So, is it possible to increase
it by factor and thereby over relax it than what we can do and that is the essence of
this. So, suppose we write x k plus 1 equal to x
k plus omega times delta x k, where delta x k is what we get using our conversions scheme.
So, if you put omega equal to 1 then, we have no over relaxation. So, we are saying that
x k plus 1 is the previous value and whatever your Jacobi method or Gauss-Seidel method
praise this delta x k to be. We are saying that it is only asymptotically approaching
and it can be higher. So, can it be higher by a factor omega? If you say omega is greater
than 1 then, you are over relaxing at every step so it becomes a successive over relaxation
method. So, this relaxation methods are by which you
can find delta x k and you are saying that is more than that by certain amount omega
when omega is greater than 1, for typical problems where you apply Gauss-Seidel method
all though not in all cases, in the case where a is diagonally dominant is good way of saying
it. In such a case you can show that omega can be less than 1 and this modification factor
will not affect it is property of whether it is converging or diverging.
So, the method will remain convergent as long as omega is between 0 and 2. If it is less
than 0 or greater than 2, then it may diverge. So that means, that we can take any value
between these and typically when omega is less than 1, you say it is under relaxation.
When omega is greater than 1 you say it is over relaxation. So when omega is greater
than 1 you are taking a step size which is greater than what should be taken as per the
Jacobi method or the Gauss-Seidel method. When omega is less than 1 you are saying that
it should be less than what it should be taken what is estimated.
If you are trying to speed up the rate of convergence, speed up the number of steps
that are required to get your objective of getting close enough to the true value, then
you would like an over relaxation omega to be greater than 1. But if you suspect, for
example this delta x k estimate itself is subject to errors because you are solving
a set of equations and you are making assumptions because of coupling or non-linearity or making
assumptions about what is the delta x k? That is when you are trying to solve non-linear
algebraic equation or a set of equations, then you could say for the sake of safety,
for the sake of convergence you could say that i will under relax i would not take the
full value that is predicted and only take part of the value.
So, then you could say that i will take omega to be between 0 and 1, you will be under relaxing,
but here are a, x equal to b is a linear equation linear algebraic equations and we do not have
to worry about under relaxation. We can definitely use over relaxation and when a is diagonal
dominant we can take omega to be anything between 1 and 2 and we can do that over relaxation.
So, once you put a value to omega here, then given that you have delta x k determined by
Jacobi method or whatever method. So, if you use Gauss-Seidel method, this is
equivalent to saying that x k equal to m x x k plus 1 equal to m x k minus 1 plus n in
that m is d by omega minus e and n is 1 minus omega by omega times d plus f. Where, d is
the diagonal terms of a and e is the those terms which are below the diagonal and f for
the terms which are above the diagonal. So, with this notation for a given a in a x equal
to b you can directly find m and n as per the Gauss-Seidel method with this extra factor
omega which is added and we take for over relaxation omega to be between 1 and 2 .
So, this is will give us an iteration scheme like this. So, now you have x k plus1equal
to d by omega minus e inverse of this whole thing here and what we would like to say is
that in this case if you put this as x k plus 1 equal to p x k plus q, this p here is different
from the p that you have with g s with the Gauss-Seidel method. What it means is that,
it is convergence rate is affected by the change that you are making and the change
that you are making is omega here. If you take omega to be 1 then this thing
will go to 0 and you will have d by 1 is d and you will get back to Gauss-Seidel method.
If omega is greater than 1, then you have extra influence of the diagonal elements here
and these things here. So, that will mean that you have and iteration matrix, which
is different from what you have it is the Gauss-Seidel method and therefore, it will
have different convergence behavior. Now all though we are putting it like this
in terms of this inverses as before, as in the case of Jacobi method and Gauss-Seidel
method, we do not do matrix inversion directly. This whole thing can be written just like
we wrote in the Gauss-Seidel method except that, we put this a11 x 1 k plus 1 here and
then we have this a 1 1 x 1 k minus omega times this whole thing and what is the omega
times a 1 1 x 1 k plus a 1 2 x 2 play that is all the elements at the previous time step
minus b 1and this whole thing is minus b 1. So, this is b 1 times omega minus this value
and plus this value here. And similarly for x 2 k it is a 2 2 x 2 k
minus omega times this whole thing. Since we know x 1 as the updated value we make use
of the updated value, as for as evaluating as per the successive over relaxation is concerned,
we are adding 1 multiplication here. Omega times this and then this is we can depending
on how it is computed. So, it is either already computed or it is there and this gives a subtraction.
So, and then this whole thing is divided by a 1 1.
So, I can see as a result of this may be 1 or 2 more multiplications per equation and
since you have equation here is like that for us pass matrix you have only certain number
of non 0 components here. The total number of arithmetic operations to go from step k
to step k plus 1, is not very different from that for Jacobi or Gauss-Seidel method and
it varies as may be six n or seven n and like that compare to 5 n for the case are Gauss-Seidel
method or Jacobi method. So, despite putting this extra factor here,
you are only increasing the number of arithmetic operations by n in this. So, it becomes and
instead of factor of 5 you now have factor of 6 or 7 depending on how you compute it
here. That makes it readily implementable as easily implementable, as the Gauss-Seidel
method or Jacobi method. So, in a way it is a small extension to the Gauss-Seidel method.
Just as Gauss-Seidel method can be considered as a small extension to the Jacobi method.
So, this is all in the same family with the same kind of approach solution. Now what does
it did to the convergence rate? We have said that as long as this omega is between 0 and
2, whether it diverges; converges is not changed. Convergence is assured, but at what rate?
Now the rate of convergence is depends as before on the spectral radius of the new iteration
matrix. So, that is involving all this omegas and
it cannot be computed for the general case. For any case of either that may be resulting
from the discretization like this, but for the specific case of a symmetric and positive
definite. So, then optimum value of omega, so the idea
is that this omega varies between 1 and 2 for the case of over relaxation, but the rate
of convergence does not vary monotonically between 1 and 2, it typically goes through
a minimum or the number of iterations needed to reduce a residual by factor of 10. That
value which was saying is asymptotically convergence rate. It reduces as omega is increased it
goes to a minimum and then it is starts again coming back up. So that means, there is an
optimum value of omega at which for a given a matrix and for a given SOR method.
There is in optimum value for which you will have the least number of iterations steps
need to reduce the residual by factor of 10. So, it will have the highest convergence rate.
What the optimum value is not known apriori, but for the specific case of when a is symmetric
in positive definite, then the optimum value for Gauss-Seidel method of SOR where delta
x k is determined as per the Gauss-Seidel method is given by this expression here. So,
this is 2 divided by 1 minus 1 minus rho square root of 1 minus rho square and what is this
rho? it is rho is a spectral density of the optimum value of optimum iteration matrix
think. So, where the Jacobi s o r is optimal, when
this is known when the optimal value of omega is known, then the spectral radius is defined
as omega optimal minus 1. So, this is a spectral radius. For the optimal value of Gauss-Seidel
SOR and this it depends on the spectral radius of this. For the specific case of Laplace equation with Dirichlet boundary conditions, you know
all this things. We know the Eigen values and we know everything about it and the spectral
radius of the Jacobi iteration matrix is cosine pi by m therefore, the optimum value of omega,
the over relaxation parameter for the Gauss-Seidel SOR scheme is given by 2 by 1 plus sine pi
m which is roughly equal to 2 times 1 minus pi by m and therefore, the spectral radius
is given by this minus one. So, that gives us 1 minus 2 pi by m.
Now, how does this is compared with the previous values? For the previous value for the g s
scheme we had rho g s as 1 minus pi square by m square. What do we have here? We have 1 minus 2 pi
by m and what is significant is that in the case of the spectral radius for Gauss-Seidel
it is pi square by m square and here it is m and what is m? m is the number of divisions
in the x direction and n the number of equations is the number of divisions in the x direction
times the number of divisions in the y direction because that gives you the total number of
grid points. So, in the case where the number divisions
in the x direction and y direction are the same, this m is equal to square root of n,
where as in the case of Gauss-Seidel method without SOR m is you had m square here so,
it is proportion to m. Now what does that mean? That for large m this 2 pi by m is smaller
than pi square by m square so that means that, and the spectral radius here is 1 minus this
value. So, this value is higher than the spectral radius is smaller. When spectral radius is
smaller the rate of convergence is faster. So, for large values of m therefore we can
expect the successive over relaxation method of the Gauss-Seidel think to be faster than,
what it is for the simple Gauss-Seidel method. So, for large value of m and in the asymptotic
limit, the number of iterations require to reduce error or residual by in order magnitude
varies as square root of n, where n is the number of equations or number of grid points
and the total number of arithmetic operations required to reduce error by a factor of 10
varies as, n raise power 3 by 2 or enter power 1.5 and this value varies as n square for
Gauss-Seidel method so, that means that the arithmetic number of multiplications or divisions,
is reduce from constant times n square to constant times n to power 1.5 or square root
of n cubed. How significant is it? If you take million
grid points then n square is million square. So, 10 to the power 12 and n to the power
1.5 is million times square root of million. So, that is 10 to the power 3. So, that is
10 to the power 9. So, that means, that instead of taking 10 to the power 12 number of multiplications
you are taking 10 power 9 number of multiplications. So, that is multiplicative factor of 1000.
So that means that, s o r method is almost 1000 times faster than the non SOR Gauss Seidel
method. So, that is a kind of computation advantage we can get, if we were using the
optimal value of omega. Optimal value of omega is known only for set in special cases.
If you want to find out for the true case then you have to do much more and you have
to known the all the Eigen values, and the determination of Eigen values for a general
matrix itself will take n cube number of mathematical operations so that means that, we cannot look
at we cannot say that let us take the matrix and let us find the Eigen values and then
try to choose try to find the best optimal value all that is not possible because Eigen
value determination is take itself will take n cube number of operations. So, that is why
only for certain class of problems under certain cases we know it is a convergence rate and
we know that it can converge very fast by significant margin compares to the Gauss-Seidel
method. So, in a general case the optimum value of
omega changes with a, with the matrix. So,1 may need to estimate numerically. So, that
is you do for 50 equations and you see by how much by what factor the residual has decreased
for a given value of omega. For example, you take omega to be 1.5 and now you do another
50 with omega equal to1.6 and see whether that residual reduction factor has increased
or decreased. If their residual reduction factor rate is decreased, that is; if it has
taken the amount of reduction is more than what you had with 1.5 so then that means that
may be you can go to 1.7 and then you keep on going like that until you hit the reverse
trend. So, that is now you take 1.8 and you find
that it has not reduced in the 50 iterations by as much as it has reduce with value 1.7.
Then the optimal values between 1.7 and 1.8, now you take you do another 50 iterations
with 1.75 and then try to locate that kind of thing and then you can finally, get an
estimate if not the exact value you can get an estimate value of this and then use that
value to continue to drive the residual down to you are decide value.
So, that is how you can make use of this. But all this is in a way theoretical because
it takes certain number of iterations before you can get into the asymptotic convergence
limit and it is under the asymptotic convergence limit, will you have good success with the
determination of the optimal value. So, the overall margin may not be like 1000
fold it may be like 10 fold it may be like 5 fold, but still it is definitely worth it.
So, SOR method is a simple extension of the Gauss-Seidel method and for the class of Laplaces
type of equations, Parson type of equations, diagonal dominant conditions it is it can
be use very quickly very effectively to improve the rate of convergence.
In the next lecture we look at other methods having different kinds of philosophies for
increasing rate of convergence over these basic methods.
Thank you.

Jerry Heath

One Comment

Leave a Reply

Your email address will not be published. Required fields are marked *