Introduction to Iterative Methods -- Lecture Material

  
  2-D Laplace Difference Equation (LDE)
        
                        U(m 1,n)+ U(m,n 1) + 4U(m,n) U(m+1,n) U(m,n+1) = O
                
                                      m = 1,mmax ; n = 1,nmax
                                 
        Boundary Conditions:
        
                           U(O,n) = B1(n) ; U(mmaxs1,n) = B3(n) , n = 1,nmax

                           U(m,O) = B2(m) ; U(m,nmaxs1) = B4(m) , m = 1,mmax

                
        This represents a system of mmax*nmax linear equations. Such systems

        are frequently encountered in many computational disciplines.

        
        There are two basic approaches to solving such linear systems:
        
          1. Direct Methods
          
          2. Iterative Methods


Direct Methods (Gaussian Elimination and its variants)

1. Require we formulate the problem in the linear algebra format

                             AV = b
                             
2. Using a forward/back substitution (possibly with pivoting and scaling)

   the exact solution V ( up to roundoff error) is obtained in a finite number

   of arithmetic operations

   
Main Advantage:

1. The method is "direct" and straightforward -- it is clear when you are fin-
ished with the solution process.

Main Disadvantage:

1. For large numbers of linear equations, memory requirements may grow
out of bound.

2. Difficult to implement efficiently on vector and parallel architectures




Iterative Methods

1. Do not require we formulate the problem in the linear algebra format

                              AV = b
                              
2. A sequence of successive approximations to the exact solution of the lin-
ear system is generated.

Main Advantages:

1. Uses less computer memory than direct methods

2. Easy to program for vector and parallel architectures.

Main Disadvantage:

1. The method is not "direct" -- it is not clear when you are finished with the
solution process.

Introduction to Iterative Methods:

Recall the 2-D Laplace Difference Equation (LDE):


      - U(m-1 ,n) - U(m,n-1) + 4U(m,n) - U(m+1,n) - U(m,n+1) = O
      
                    m = 1,mmax ; n = 1,nmax
                    
We can reformulate this as

        U(m,n) = 25*[U(m-l ,n) + U(m,n-1 ) + U(m+1 ,n) + U(m,n +1 )]
        
and then, starting from some initial guess of the solution, generate succes-
sive approximations according to:

 U^(k+1)(m ,n) = .25*[U^k(m,n-1 ) + U^k(m+1 ,n) + U^k(m ,n+1 )]
       
                         m = 1,mmax ; n = 1,nmax

Note:  U^(k+1) is to be read as "U superscript (k+1)" and U^k is "U superscript k"
                         
This is called the (point) Jacobi Iterative Method.

How would you code this method?? Would your code vectorize??

Code Fragments for Gauss-Seidel Method

    --------     Initialize and set boundary conditions for U --------
            
        do 100 k = 1 ,kmax

        do 10 n = 1,nmax

        do 20 m = 1 ,mmax

        U(m,n) = .25*[U(m-1,n) + U(m,n-1) + U(m+1,n) + U(m,n+1)]

        
  20 continue
  
  1 0 continue
  
                resmax = 0.0

                do 1 n = 1,nmax

                do 2 m = 1,mmax

                res = 4*U(m,n) - [U(m-1,n) + U(m,n-1) + U(m+1,n) + U(m,n+1)]

                resmax = amaxl(abs(res),resmax)

        
         2     continue
  
        1       continue

                if(resmax .LT. tol) go to 300 ! tol is a preset closure criterion

         
  100 continue


  




            Note: On a problem on the unit square with
       
                     a 50 x 50 interior grid

                     tol set to .0005

                     boundary condtions set to x^2 -y^2

                     
       the number of iterations required for convergence is
       
           Jacobi               829

           Gauss Seidel     420

           Optimal SOR*  103

           
       * The optimal value of the relaxation parameter for this problem is
       
                 omega = 1/(1 + sin(pi/51 ))