Unit 8: Numerical Methods

Chapter 39: Systems of Equations

Section 39.6: relaxation and SOR

Copyright

Copyright * 2001 by Addison Wesley Longman, Inc.

All rights reserved. No part of this publication may be reproduced, stored in a retrieval system, or transmitted, in any form or by any means, electronic, mechanical, photocopying, recording, or otherwise, without the prior written permission of the publisher. Printed in the United States of America.

Initializations

> restart;

> with(linalg):
with(plots):
read(`C:/Program Files/Maple 6.01/pvac.txt`):

Warning, the protected names norm and trace have been redefined and unprotected

Warning, the name changecoords has been redefined

>

Relaxation - the Concept

The idea of relaxation has been attributed to the British engineer, Richard Southwell, [1], although some claim that Gauss was familiar with the concept. Conceived in the "by-hand" era before computers, the idea underlies the more modern method of successive overrelaxation (SOR). Our exposition will be in terms of the linear system A x = y for which the matrix A is the strictly diagonally dominant

> A := matrix(3,3,[-12,-1,1,1,13,7,3,4,8]);

>

and for which the vector y is

> y := vector([-11,48,35]);

>

Then, the solution is

> xe := linsolve(A,y);

>

the system variables are

> vars := [seq(x[k],k=1..3)];

>

and the equations of the system are

> for k from 1 to 3 do
q||k := evalm(A &* vars)[k] = y[k];
od;

>

The residual vector for this system would be R = A x y , which we make into a function in Maple via

> R := map(unapply,evalm(A&*vars - y),op(vars));

>

Clearly, R would be zero if x = x . Indeed, we find

> R(op(convert(xe,list)));

>

If we start an iteration process at x = 0 , the residual vector would initially be R = y , in which the largest absolute value is 48, occurring in the second equation. In fact, we see that

> R(0,0,0);

>

If we could make the largest residual zero, that is, if we could relax the largest residual to zero , we should be led to a more accurate solution. Our strategy, then, is to keep = 0, and adjust to make the second component of R become zero. We would be aided in this task if we solved each equation for the variable with the coefficient of greatest magnitude in that equation, obtaining

which we implement as functions in Maple via

> for k from 1 to 3 do
X[k] := unapply(solve(q||k,x[k]),op(remove(has,vars,k)));
od;

>

So far, the astute reader will recognize the fixed-point iteration formula x = g ( x ). However, the actual iteration we do will be different.

The value of which drives the second component of R to zero is

> X[2](0,0);

>

If is set equal to this value while keeping = 0, then the next member of the sequence of iterates is

x = (0, , 0) =

which we represent in Maple as the list

> z := [0,X[2](0,0),0];

>

Since the vector x has changed, the residual vector R will change to

> `R` = evalf(R(op(z)));

>

Now, the third component of the residual vector has the largest magnitude, so must be varied to bring this component of the residual to zero. The appropriate value of is computed by

> new := X[3](z[1],z[2]);

>

so x , the next iterate, is given by

> z := evalf([z[1],z[2],new]);

>

The residual generated by x is then

> `R` = R(op(z));

>

The second component of R has the largest magnitude, so is adjusted to bring this component of the residual to zero. This is done by giving the value

> new := X[2](z[1],z[3]);

>

so x , the next iterate is given by

> z := [z[1],new,z[3]];

>

The residual generated by x is then

> `R` = R(op(z));

>

The first component of the residual now has the largest magnitude, so must be adjusted to bring this component of R to zero. This is done by assigning the value

> new := X[1](z[2],z[3]);

>

so x , the next iterate, is given by

> z := [new,z[2],z[3]];

>

The residual generated by x is then

> `R` = R(op(z));

>

The component of R with the largest absolute value is the third, so must be adjusted to bring this component of the residual to zero. This is done by assigning the value

> new := X[3](z[1],z[2]);

>

so x , the next iterate, is given by

> z := [z[1],z[2],new];

>

The residual generated by x is then

> `R` = R(op(z));

>

The second component of R now has the largest magnitude, so must be adjusted to bring this component of the residual vector to zero. This is done by assigning the value

> new := X[2](z[1],z[3]);

>

so that x , the next iterate, is given by

> z := [z[1],new,z[3]];

>

The residual generated by x is then

> `R` = R(op(z));

>

in which the third component has the largest magnitude. We therefore drive this third component of R to zero by assigning to the variable the value

> new := X[3](z[1],z[2]);

>

giving us as the next iterate the vector x , represented in Maple as

> z := [z[1],z[2],new];

>

The residual vector generated by x is then

> `R` = R(op(z));

>

The second component of R now has the largest magnitude, so it is driven to zero by assigning to the value

> new := X[2](z[1],z[3]);

>

This gives x , the next iterate, as

> z:=[z[1],new,z[3]];

>

and the iteration continues until it converges to the exact solution

> print(xe);

>

Although it can become tedious, selecting the component with the largest magnitude in the residual vector is something which can be done "by-hand" for a system of modest size. Clearly, having to pick, repeatedly, the largest number from a list of 10,000 would be highly susceptible to error. Moreover, having a computer make the many comparisons required to select the largest magnitude is not efficient, and this version of the relaxation algorithm is not a viable method for solving large problems on a computer.

However, the underlying concept does motivate the method of successive overrelaxation, the SOR technique, which is capable of being implemented on a computer.

>

Successive Overrelaxation

Southwell discovered that often, when driving a component of the residual to zero, the other components would move in the same direction. For example, the residual generated by x is

R =

The third component is driven to zero by increasing it from to 0. The residual generated by x is

R =

Driving the third component of the residual upwards from to 0 also drove the first component upwards from 7.3 to 9.8, and the second component upwards from 0 to 17.7.

This led Southwell to "overcompensate" or overrelax , pushing the residual, not to zero, but beyond zero. All that is needed is a more efficient way of implementing the calculations on a computer. Formulating the system A x = y as if for performing Gauss-Seidel iteration is the basis for this efficiency.

Assuming that the th equation has been solved for , we write the linear system in the form

as in the Jacobi method. Then, recognizing that the expression in parentheses is actually the th component of the residual vector R , we write

R

By introducing , the relaxation factor, this equation is modified to

R

Finally, write the general equation in the form

which is now consonant with the Gauss-Seidel method, except for the introduction of , the relaxation factor. Clearly, if , the SOR technique is precisely the Gauss-Seidel method. If A is a symmetric positive definite matrix and , then the SOR technique converges from any initial vector. (See [2] for a proof.)

>

Implementation

Implementing the SOR method for small systems is more conveniently done in matrix form. To accomplish this, write the linear system A x = y in the form

y A x = 0

and multiply though by , giving

0 = ( y A x)

Write A = L + D + U , where L is zero on and above the main diagonal, U is zero on and below the main diagonal, and D is zero off the main diagonal. In addition, add D x to each side of the equation, thereby obtaining

D x = D x ( L + D + U ) x + y

Bring the term L x to the left, giving

( D + L ) x = [ D U ] x + y

Let the vector x on the left side be the updated one, namely, x , while keeping the x on the right side as the "old" one, namely, x . This gives

( D + L ) x = [ D U ] x + y

Solve for x to obtain

x = ( D + L) [ D U ] x + ( D + L ) y

which has the fixed-point form

x = B x + C = g ( x )

where

B = ( D + L) [ D U ]

and

C = ( D + L ) y

>

Example 39.12

Consider the linear system A x = y for which A is the symmetric, positive definite matrix



> A := matrix(3,3,[8,-5,-7,-5,14,1,-7,1,8]);

>

and y is the vector

> y := vector([-23,26,19]);

>

We see that A is symmetric. That it is positive definite also, we see from

> definite(A,positive_def);

>

or from the actual eigenvalues

> eigenvals(evalf(op(A)));

>

which are clearly all positive.

the exact solution is

> xe := linsolve(A,y);

>

Decompose A as L + D + U , where L , D , and U are given by

> L := matrix(3,3, (i,j) -> if i>j then A[i,j] else 0 fi);
U := matrix(3,3, (i,j) -> if i<j then A[i,j] else 0 fi);
d := diag(seq(A[k,k],k=1..3));

>

Notice that we have used d in place of D , which in Maple is reserved for the differentiation operator.

The matrix B = ( D + L ) [ D U ] is then

> B := map(simplify,evalm(inverse(d+omega*L)&*((1-omega)*d-omega*U)));

>

and the vector C = ( D + L) y is then

> C := evalm(omega*inverse(d+omega*L)&*y);

>

The eigenvalues of B determine the rate of convergence of the SOR iteration. In terms of , they are

> q := eigenvals(B);

>

Since we wish to pick so that the largest magnitude of the eigenvalues is as small as possible, we begin with Figure 39.5, a graph of the three eigenvalues as a function of .

> p1 := plot(map(abs,[q]),omega=0..2,color=[black,red,green], thickness=[3,1,1], labels=[``,``]):
p2 := textplot([1.75,-.07,`w`], font=[SYMBOL,12]):
display([p1,p2],xtickmarks=5, ytickmarks=5);

>

The optimal value of appears to be approximately 1.55, determined from the graph. For the matrix B becomes

> B1 := subs(omega=1.55,op(B));

>

and the vector C becomes

> C1 := subs(omega=1.55,op(C));

>

Starting with the initial vector , we obtain the first ten iterates and the norms of the differences between the iterate and the exact solution.

> x := evalf(vector([10,10,10]));
for k from 1 to 10 do
x := evalm(B1&*x+C1);
print(k,x,norm(x-xe));
od:

>

To verify that the iteration converges more slowly for other values of , we try and obtain

> B2 := subs(omega=1.75,op(B));

>

for the matrix B , and

> C2 := subs(omega=1.75,op(C));

>

for the vector C . Again iterating from the initial vector , we find

> x := evalf(vector([10,10,10]));
for k from 1 to 10 do
x := evalm(B2&*x+C2);
print(k,x,norm(x-xe));
od:

>

With , the tenth iterate is farther away from the exact solution than with . We conclude that our choice of based on the graph of the magnitudes of the eigenvalues was appropriate.

>

Equivalence of Forms

When , the SOR method reduces to the Gauss-Seidel method, an equivalence we will now verify for the system in Example 39.12 whose Gauss-Seldel equations are

When , the matrix B becomes

> BB := subs(omega=1,op(B));

>

and the vector C becomes

> CC := subs(omega=1,op(C));

>

Consequently, the SOR equations in the fixed-point iteration x = B x + C are

and we want to show these equations are indeed the Gauss-Seidel equations. Clearly, the first equation in each formulation is the same. In the second Gauss-Seidel equation, replace from the first, obtaining

the second equation in the SOR formulation. In the third Gauss-Seidel equation, replace from the second equation. This introduces another . Now replace all from the first equation, and the result will be

the third equation of the SOR formulation. Thus, we have verified that the matrix form of the SOR iteration actually uses the updated values x in the same way that Gauss-Seidel iteration does..

To execute these substitutions in Maple, write the Gauss-Seidel equations as

> unassign('x','X');
for k from 1 to 3 do
X[k] := solve(evalm(A &* vars-y)[k], x[k]);
od;

>

Substitute for in the second equation, using its value from the first equation, thereby obtaining

> X2 := subs(x[1]=X[1],X[2]);

>

This is the second equation in the SOR formulation.

Finally, in the third equation of the Gauss-Seidel formulation the variables and both appear. The variable must be updated with the most recent version of this variable. The variable must likewise be replaced with its updated version, but that update contains the old value of . Thus, the replacement of must be done after the replacement of , lest the in be updated twice. With this due care, we get

> subs([x[1]=X[1], x[2]=X2], X[3]);

>

on the righthand side of the third equation, and we have verified that the matrix form of the SOR iteration actually uses the updated values x in the same way that Gauss-Seidel iteration does.

>

References

[1] Richard V. Southwell, Relaxation Methods in Engineering Science, Oxford University Press, 1940.

[2] James M. Ortega, Numerical Analysis, A Second Course, Academic Press, 1992.

>