Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Solving a linear system with Gauss-Seidel's method

  1. Jun 27, 2011 #1

    fluidistic

    User Avatar
    Gold Member

    Edit: Problem solved. Look at the bottom of this post.
    Ok so far I've made a working program that can solve for x in Ax=b using Jacobi method.
    I've tried to use Gauss-Seidel method too, but I failed in doing so.
    I took as an example, the following A matrix:
    Code (Text):
    4 -1 0 -1 0 0
    -1 4 -1 0 -1 0
    0 -1 4 0 0 -1
    -1 0 0 4 -1 0
    0 -1 0 -1 4 -1
    0 0 -1 0 -1 4
     
    with the following b vector:
    Code (Text):
    1
    1
    1
    1
    1
    1
    I use the initial "guess" or vector:
    Code (Text):
    1
    1
    1
    1
    1
    1
    .
    As a result, I should get a value for the vector x close to
    Code (Text):
    0.57024383544921875    
      0.71393871307373047    
      0.57024383544921875    
      0.57118320465087891    
      0.71261024475097656    
      0.57118320465087891  
    but I get the values
    Code (Text):
    0.27777910232543945    
      0.62212115526199341    
      0.54369938373565674    
      0.48600935935974121    
      0.66520172357559204    
      0.55222527682781219    
    Here comes my Jacobi and then Gauss-Seidel codes inside the same program:
    Code (Text):
    program hope
    implicit none

    integer :: p, n, n_max
    real :: tol
    98 write(*,*)"In order to solve the linear system, do you want to use Jacobi(1)'s method or Gauss-Seidel(2)'s one?"
    read(*,*)p
    write(*,*)"Enter the order n of the matrix A, a tolerance and a maximum number of iterations"
    read(*,*)n, tol, n_max
    if (p==1) then
    call jacobi(n_max,n,tol)
    else if (p==2) then
    call gauss(n_max,n,tol)
    else
    write(*,*)"Please choose a valid option"
    goto 98
    end if
    end program


    subroutine jacobi (n_max,n,tol)
    implicit none
    integer                             :: i,j,k
    real                                :: sup_sum, low_sum, norm
    integer, intent(in)                 :: n_max, n
    real, dimension(n)                  :: b, u,x
    real, intent(in)                    :: tol
    real, dimension(n,n)                :: A  



    open(unit=10,file="A_matrix.dat")
    open(unit=15,file="b_vector.dat")
    open(unit=20,file="initial_vector.dat")
    do i=1,n
      read(10,*) (A(i,j),j=1,n)
      read(15,*) b(i)
      read(20,*) x(i)
    end do
    close(10)
    close(15)
    close(20)
    k=1
    norm=100
    do while ( (k<=n_max).and.(tol<norm) )

      sup_sum = b(1)
      do j=2,n
        sup_sum = sup_sum - A(1,j)*x(j)
      end do
      u(1)=sup_sum/a(1,1)
     
      do i=2,n
        low_sum=dot_product( A(i,1:i-1), x(1:i-1))
        sup_sum=dot_product( A(i,i+1:n), x(i+1:n))
        u(i)=(b(i)-low_sum-sup_sum)/A(i,i)  
      end do
     
     
      norm=0
      do i=1,n
      norm=norm+abs(x(i)-u(i))
      end do
     
      do i=1,n
        x(i)=u(i)
      end do
      write(*,*) k, norm
      k=k+1
    end do

    open(unit=25,file="final_vectorjacobi.dat")
    do i=1,n
      write(25,*) x(i)
    end do

    close(25)


    end subroutine


    subroutine gauss (n_max,n,tol)
    implicit none
    integer                             :: i, j, k
    real                                :: low_sum, sup_sum, norm
    real, dimension(n,n)                :: A
    real, dimension(n)                  :: b, x, oldx
    integer, intent(in)                 :: n_max
    real, intent(in)                    :: tol
    integer, intent(in)                 :: n
     


    open(unit=10,file="A_matrix.dat")
    open(unit=15,file="b_vector.dat")
    open(unit=20,file="initial_vector.dat")

    do i=1,n
      read(10,*) (A(i,j),j=1,n)
      read(15,*) b(i)
      read(20,*) x(i)
    end do
    close(10)
    close(15)
    close(20)

    k=1
    norm=100
    x=0
    do while ( (k<=n_max).and.(tol<norm) )
    oldx=x

      do j=2,n
        sup_sum = sup_sum - A(1,j)*x(j)
      end do
      x(1)=sup_sum/a(1,1)

      do i=2,n
        low_sum=dot_product( A(i,1:i-1), x(1:i-1))
        sup_sum=dot_product( A(i,i+1:n), x(i+1:n))
        x(i)=(b(i)-low_sum-sup_sum)/A(i,i)  
      end do
     
      norm=0

     
      norm=sum(abs(x-oldx))
        write(*,*) k, norm
       
    k=k+1
    end do
    open(unit=25,file="final_vectorgauss.dat")
    do i=1,n
      write(25,*) x(i)
    end do
    close(25)
    end subroutine
    I must say I'm a bit confused in Gauss-Seidel subroutine when I treat x and oldx as arrays rather than real numbers. I know I do get wrong values for x(i) for Gauss-Seidel method... but I don't see the error.
    I'd appreciate any help to search the error(s). Thanking you in advance.


    EDIT: Sorry guys! I think I've found out the error by a pure lucky guess. I changed the 2 starting values for the loops from 1 to n instead of from 2 to n. Well I won't erase this thread, so that everyone can use this code :)
     
  2. jcsd
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook

Can you offer guidance or do you also need help?
Draft saved Draft deleted



Similar Discussions: Solving a linear system with Gauss-Seidel's method
Loading...