How to Solve a System of Linear Equations using the Gauss-Seidal Method

The Gauss-Seidal Method is another iterative method used to solve a system of linear equations. It is also known as the method of successive displacement.It converges only if the matrix is strictly diagonally dominant or symmetric and positive definite. Let us see how to use this method to solve linear equations.

• We will need to give the number of variables or the size of the coefficient matrix. The column order will be m=n+2. This is because we will add two more columns, one that will be part of the augmented matrix(the R.H.S or constants in the equations) and another which will store the initial guesses for the solutions x1 through xn.
• Next we set maxcount, the maximum number iterations, at the terminal. Then we read the input file and display it on the screen.
• We call a subroutine to solve the equation. The subroutine first tests the coefficient matrix for strict diagonal dependence. And , gives a go ahead only if it satisfies this condition. Else, the method fails. To test for diagonal dependence we simply sum the absolute values of all elements other than those on the diagonal for every row. We then test if the absolute value of diagonal elements in each row is greater than the sum of absolute values of all other elements in that row.
• We set an array b(:)=a(:,n+1), which is the R.H.S of the actual equations. And we set the array x=0, our initial guesses.
• In a DO loop, we set oldx=0 and start changing values of x(i) using two more DO loops. At the end of these loops , we check  diff=abs(x-oldx) and error=sum(diff). If error<eps, we return and display results in the main program. If not, we continue running the two DO loops till this is attained.
• Sometimes, the method does not converge. In this case we display "maximum number of iterations exceeded" if pcount>maxcount.

Let us understand the process better with an illustration.

Gauss-Seidal Method for a System of 3 Linear Equations

Now for the Fortran Program,


    module seidal
        real, allocatable, dimension (:) :: a(:,:),x(:),b(:)
    end module seidal

!To find solutions for a system of linear equations using the Gauss Seidal Method
program seidal_main
    use seidal
    implicit none
    integer :: m,i,j,flag,n,maxcount
   
    write(*,*) "Give the row order of matrix"
    read (*,*) n
    write(*,*) "The row order of the matrix is",n

    m=n+2
    write(*,*) "The column order of the matrix is ",m

      write(*,*) !To read the maximum number of iterations
    write(*,*) "Give the maximum count of iterations"
    read (*,*) maxcount
    write (*,*) "The maximum number of iterations is", maxcount
   
    allocate (a(n,m),x(n),b(n))
    write(*,*) "The augmented matrix is given below along with initial estimates:"
   
    open (unit=1,file="gs111.dat")
   
    do i=1,n
          read (1,*) (a(i,j),j=1,m)
      end do
   
    do i=1,n
          write (*,*) (a(i,j),j=1,m)
      end do

    call linear (m,flag,maxcount)
       if (flag==0) then
          write(*,*) "The coefficient matrix does not possess strict diagonal dominance. Hence, the Gause Seidal method fails."
          else
        write(*,*) "The matrix containing x values i.e. x1...xn"
        do i=1,n
            write(*,*) x(i)
        end do
          write(*,*) "The x values obtained satisfy the equations"
    end if

end program seidal_main
         
    subroutine linear (m,flag,maxcount)
        use seidal
        implicit none
       
        integer:: i,j,flag, pcount,m,maxcount,n
        real:: error,eps,sumrow,absD,sumND,sumax
         real, allocatable::oldx(:), diff(:)
       
        eps=10.0e-15
        n=m-2
        allocate (oldx(n),diff(n))
        do i=1,n !To check the validity of the method (strict diagonal dominance)
            sumrow=sum(abs(a(i,1:n)))
            absD=abs(a(i,i))
            sumND=sumrow-absD
            if  (absD.le.sumND) then
                  flag=0
                  return
                  else
                flag=1
            end if
        end do
       
        pcount=0
        b(:)=a(:,n+1)
        x=0
        do !To begin the iteration
            oldx=x
           
            do i=1,n
                sumax=0.0
               
                do j=1,n
                    if (j<i) then
                        sumax=sumax+a(i,j)*x(j)
                    end if
                    if (j>i) then
                        sumax=sumax+a(i,j)*oldx(j)
                    end if
                end do
               
                x(i)=(b(i)-sumax)/a(i,i)
            end do
           
            diff=abs(x-oldx) !To find the sum of errors
            error=sum(diff)
            if (error<eps) exit
 
                   pcount=pcount+1
                if (pcount>maxcount) then
                    write(*,*)  "No. of iterations exceeded the maximum no. of iterations"
                       stop
                end if
        end do
    end subroutine linear
                         

Sample Output for Gauss Seidal Method for a System of 3 Linear Equations