Gauss-Seidel method

หลังจาก Jacobi ได้นำเสนอวิธการแก้ระบบสมการ มาระยะหนึ่ง ก็มีคนนำมาพัฒนาต่อ แต่หลักการจะใกล้เคียงกัน

สมมติเรามีระบบสมการเชิงเส้น  3 สมการดังนี้

จัดรูปสมการใหม่ตามหลักการ Jacobi method ได้ดังนี้

เมื่อต้องการหาค่า x1 ก็ให้ลองใส่ค่า x2 , x3 ในสมการ แล้วนำค่า x1 ใหม่ดังกล่าวไปใช้ในการคำนวณค่า x2 เลย และก็นำ x1,x2 ใหม่ที่เพิ่งคำนวณได้ไปใช้ในการหา x3 ได้ทันที จนถึง xn จะเห็นว่า Gauss-Seidel method มีจุดต่างกับ Jacobi คือ xi ใดๆ ที่เพิ่งหามาได้จะถูกนำไปแทนค่า ในการหาค่า ของตัวแปร Unknown อื่นๆทันที ในรอบนั้นเลย และทำเช่นนี้ไปเรื่อยๆจนค่า x1 รอบปัจจุบัน กับรอบที่ผ่านมามีค่าแตกต่างกันน้อยมากๆ หรือไม่ต่างกันเลย (ดังสมการข้างล่างนี้)  และ x2,x3 .. xn ก็เช่นเดียวกัน ถ้าเป็นเช่นนั้น ค่า x1,x2,x3 ... xn ที่ได้จากรอบล่าสุด คือคำตอบสุดท้าย

ในขณะที่ Jacobi method จะใช่ค่าเริ่มต้นใหม่ ทุกตัวในรอบต่อๆไป เท่านั้นแต่ในระหว่างรอบ จะไม่มีการนำค่าที่คำนวณได้ใหม่ มาใช้แต่อย่างใด ซึ่งจะเห็นว่าความแตกต่างดังกล่าวจะทำให้ Gauss-Seidel method จะใช้จำนวณรอบในการคำนวณน้อยกว่า Jacobi method มากเลยทีเดียว ซึ่งหากเราเขียนโปรแกรมสั่งให้คอมพิวเตอร์หาค่าให้เรา ถ้าจำนวนสมการเป็นหลัก พัน หรือ มากกว่า เราก็จะใช้เวลา หน่วยความจำของคอมพิวเตอร์น้อยกว่าการใช้ Jacobi method อย่างมากทีเดียว

ใน การแก้สมการ ผู้เขียนใช้วิธี แยกหา x1 และ xn แยกออกมาต่างหาก ส่วน x2 จนถึง xn-1 จะหาโดยใช้ loop เดียวกัน

Gauss-Seidel m-file.

function gauss_seidel(A,b,error,max)
% Use to solve the liear equation system by Gauss-Seidel's technique.
% A : NxN matrix
% b : 1xN matrix
% error : Constant number that use to check and stop computation.
% max : maximum loop that the computer will compute the result.
% General form : Ax = b
% |a11 + a12 + .. a1n | | X1 | | b1 |
% |a21 + a22 + .. a2n | | X2 | | b2 |
% . . . . = .
% . . . . .
% |an1 + an2 + .. ann | | Xn | | bn |
%
% Entry example
% A = [2,-3,1;4,-3,5;-2,1,4];
% b = [9;5;6];
% gauss_seidel(A,b,0.0001,100)
% If you type the above lines on the command windows and then press ENTER,
% you will see the X Matrix that make the equation be true.
% Wrote By :
% Chalong Srikaewsiew , Master degree of Mechatronics Engineering
% Suranaree University of Technology , Nakornratchasima, Thailand
% October 29,2005
clc
N = length(A); % Size of the matrix
% Initial value of x
  for k=1:1:N
    x0(k)= 0;
  end
% Counter and index
   count = 0;
   loop = 0;
   flag = 0;
% Start iteration
   while loop < max
       for i=1:1:N
           s=b(i);
           f=A(i,i);
           A(i,i)=0;
               for j=1:1:N
                   s=s-A(i,j)*x0(j);
               end
           s=s/f;
           A(i,i)=f;
               if abs(s-x0(i)) > error
                  x0(i)=s;
               else
                  flag = 1;
                  loop = max;
               end
      end
          count=count+1;
    end

% Report the result...
  if flag ==1
     fprintf('\t X Matrix \n');
     disp('--------------------');
     table(:,1)=x0;
     fprintf('\t %4.4f \n',table');
     disp('--------------------');
     fprintf(' Number of iterations : %d\n\n',count);
   else
     disp('No convergence, hit the maximum loop..');
     disp('You may increase your maximum loop number..');
        if max > 200
           disp('This lenear equation system may not converge');
        end
    end

 

Step to entry the linear system equation.

>> A=[2,1,1;3,5,2;2,1,4]

A =

     2     1     1

     3     5     2

     2     1     4

 >> b=[5;15;8]

 b =

     5

    15

     8

>>gauss_seidel(A,b,0.0001,50)

The result of the computation.

           X Matrix

      --------------------

             0.9996

             2.0002

             1.0002

      --------------------

      Number of iterations : 8

      >>