CS323 LECTURE NOTES - LECTURE 10

1    Systems of Linear Equations by Iterative Meth- ods

1.1    Jacobi’s Method

Suppose that we are given a system of equations:

a11x1 + a12x2 + a13x3 + ... + a1nxn    =   b1

a21x1 + a22x2 + a23x3 + ... + a2nxn    =   b2

.         .

.         .

.         .

an1x1 + an2x2 + an3x3 + ... + annxn    =   bn

we can solve the rst equation for x1  and get

b1 − (a12x2 + a13x3 + ... + a1nxn)

a11

from the second equation we can get x2

b2 − (a21x1 + a23x3 + ... + a2nxn)

22

and so on.

In general:

bi − Pjn=1;j =6i aijxj

aii

(2)

Suppose that we are given an initial guess to the solution of the system:

x,x,...,x

We can use the equations (2) to compute a new approximation to the solution of the system:

x = bi − Pjn=1;j =6i aijxj(0)

Based on this equation (Jacobi’s equation) we can write an itera- tive algorithm where we keep computing new guesses checking to see if the solution is within an ǫ of the actual solution, i.e. kA(1) −k < ǫ

1.2    Convergence

It is possible to show that Jacobi’s method will converge to a solution when kA(1) − k < ǫ i it converges when k(1) − (0) k < ǫ, which is faster to compute.

So we have the following algorithm:

Jacobi’s Algorithm

Input: A, , (0),ǫ,N

k = 0

repeat

k + +

for i = 1 to n

x = bi − Pjn=1;j =6i aijxj(0)

error=k(1)  − (0) k

(0)  = (1)

until error< ǫ or k > n

if error< ǫ

return (1)

else

throw solution_not_found_exception

Denition

A matrix A is diagonally dominant i

|aii| ≥  X  |aij |

i=1;ji

i.e.  the absolute value of a main diagonal entry is greater than or equal to the sum of the absolute values of the entire column (without including the diagonal element).

Theorem

If A is diagonally dominant then Jacobi’s method converges on matrix A

Notice that it is still possible for a matrix not to be diagonally dominant and for Jacobi’s method to converge. But if the matrix is diagonally dominant then the method converges for sure.

1.3    Gauss-Seidel’s Method

This method is a variation of Jacobi’s Method, where the newly computed values x(1)  can be used as soon as they become available.

When computing x, we know all of the new values x    j = 1 ...i − 1, so we use them, and we use x(0)  for the other values (i = j + 1 ...n).

bi − Pji −=11 aijxj(1) +Pjn=i+1aijxj(0)

aii

In general this algorithm converges faster than Jacobi.

Notice also that it is possible to reuse the  vector, in the algo- rithm, so we don’t really need to keep track of both vectors, but we need to keep track of the sum of the squares of the dierences (x − x) every time that we update the value of xi  to be able to compute the norm k(1) − (0) k.

This means that when xi  is updated we store in a variable t the previous value, and compute the dierence between the new value and the previous value squared (xi  − t)2 , and use a variable s to accumulate the summation, so that we can compute the norm as √s.