Jacobi method

เมื่อสงครามโลกครั้งที่ 2 อุบัติขึ้น สหรัฐเริ่มตระหนักถึงความแตกต่างด้านวิทยาการสมัยใหม่ระหว่างยุโรปกับอเมริกา ที่ยังห่างไกลกันและทางที่ดีที่สุดที่สหรัฐจะไล่หลังยุโรปได้ทันไม่ใช่สงคนข้ามมหาสมุทรแอตแลนติกมาเรียนที่ยุโรปอย่างที่กำลังทำอยู่ แต่เป็นการหาคนที่เก่งๆจากยุโรปไปถ่ายทอดความรู้ที่สหรัฐอเมริกา

 Jacob Jacobi เป็นนักคณิตศาสตร์ชาวเยอรมัน ที่ถือว่าโด่งดัง และมีความสามารถทางคณิตศาสตร์สูงมาก และเป็นขุมกำลังทางสมองของกองทัพเยอรมัน ในช่วงสงครามโลกครั้งที่ 2  เขาเป็นหนึ่งในเป้าหมายที่สหรัฐต้องการอยากได้ตัว เพื่อให้ย้ายไปอยู่สหรัฐอเมริกา เช่นเดียวกับ อัลเบิร์ต ไอน์ไสตน์  และภายหลังสิ้นสงครามโลกครั้งที่ 2 เขาก็ได้ไปอาศัยอยู่ในอเมริกา และก็ได้เป็นขุมกำลังสำคัญของสหรัฐอเมริกาทางด้านคณิตศาสตร์ ในโครงการสำคัญๆ ต่อมา

หนึ่งในผลงานที่สำคัญของเขา คือวิธีการการแก้ระบบสมการ ที่ปัจจุบันเราเรียกว่า Jacobi 's method ที่ใช้หาค่าตัวแปรที่ไม่รู้ค่าหลายๆตัวพร้อมกันในระบบสมการ หลักการที่สำคัญของวิธีที่เขาคิดขึ้นคือ ลองแทนค่าตัวแปรที่ไม่รู้ค่าทุกตัว แล้วหาว่าค่าที่ลองใส่นั้นทำให้สมการถูกหรือไม่ ถ้าไม่ถูกก็เปลี่ยนค่าไป ตามขั้นตอนดังนี้

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

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

เมื่อต้องการหาค่า x1 ก็ให้ลองใส่ค่า x2 , x3 ในสมการ และทำเช่นเดียวกัน นี้กับการหาค่า x2,x3  ในรอบแรกนี้

รอบต่อไปให้นำค่า x1,x2,x3 ที่ได้จากรอบแรก มาแทนค่าในสมการ และทำเช่นนี้ไปเรื่อยๆจนค่า x1 รอบปัจจุบัน กับรอบที่ผ่านมามีค่าแตกต่างกันน้อยมากๆ หรือไม่ต่างกันเลย (ดังสมการข้างล่างนี้)  และ x2,x3 ก็เช่นเดียวกัน ถ้าเป็นเช่นนั้น ค่า x1,x2,x3 ที่ได้จากรอบล่าสุด คือคำตอบสุดท้าย

เมื่อแรกค้นคิดเทคนิคนี้ขึ้นนั้นเป็นยุคที่เครื่องคอมพิวเตอร์ตัวใหญ่เท่าตึก 4 ชั้น แต่ปัจจุบันนี้ คอมพิวเตอร์พัฒนามาถึงขั้นที่เราจะสั่งทำอะไรยากๆก็ไม่ยากเกินไป ท่านอาจจะคิดว่า แค่นี้คำนวณตามแบบ Algebra ทั่วไปน่าจะง่ายกว่า สำหรับระบบสมการที่มีเพียง 3-4 สมการนั้นก็อาจจะใช่ครับ แต่ถ้าเป็นระบบสมการที่มีใช้งานจริงๆ บางครั้งประกอบไปด้วยสมการทั้งหมด เป็นหลักร้อย ถึงหลักแสนสมการ ถ้าลองมาหา Determinant , Inverse matrix  รับรองจอดแน่นอนครับ แม้แต่ใช้คอมพิวเตอร์ก็ยังยากเลยครับ ระบบสมการดังกล่าวจะเหมาะที่จะใช้เทคนิคการแทนค่าอย่างนี้มากกว่า

ในฐานะที่ผู้เขียนได้เรียนรู้ Jacobi method มา และอยากนำมาเผยแผ่ ให้ท่านผู้อ่านได้อ่าน และผู้เขียนได้เขียนโปรแกรม (M-File) สำหรับผู้ที่ใช้ โปรแกรม MATLAB เป็น ก็ลองเอาไปใช้ดูนะครับ แต่ต้องเข้าใจก่อนนะครับว่า มีบางระบบสมการเหมือนกันที่ไม่สามารถใช้ Jacobi ได้ ท่านอาจจะต้องลองค่าเริ่มต้นที่ดีด้วย เพื่อจะได้คำตอบ

Jacobi m-file.

function jacobi(A,b,error,max)

% Use to solve the linear system equation by Jacobi’s method.

% 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 |

%

%       Example

%          A = [2,-3,1;4,-3,5;-2,1,4];

%          b = [9;5;6];

%          jacobi(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.

% Write By :

%      Chalong Srikaewsiew , Master degree of Mechatronics Engineering

%      Suranaree University of Technology , Nakornratchasima, Thailand

clc

N = length(A); % Size of the matrix

for i=1:1:N

    x0(i)= 0.0; % Initial value of x

    eps(i) =1; % Initial value of condition check

end

% Counter and index

count = 1; % Counter

loop = 1; % Condition index

key = 0;

while loop < max

   for i=1:1:N % Row loop

       sum = b(i);

          for j=1:1:N % Column loop

             if i ~= j  % Row number much not the same as column number

                sum = sum - A(i,j)*x0(j);

             end

          end

       x(i)= sum /A(i,i);

       eps(i) = abs((x(i)-x0(i))/x(i));

       if eps(i) < error

           loop = max;

           key = 1;

       end

       x0(i)=x(i);

       count = count + 1;

       loop = loop+1;

   end

end

if key ==1

% Report the result...

   fprintf('\t X Matrix \n');

   disp('--------------------');

   table(:,1)=x0;

   fprintf('\t %4.4f \n',table');

   disp('--------------------');

   fprintf(' Number of iterations : %d\n',count);

else

   disp('No convergence , hit the maximum loop..');

 if max < 100

       disp('You may increase your maximum loop number..');

 else

   disp('This linear system equation may not suitable to apply Jacobi method');

   disp('Or your initial value for unknown may need to review again..');

 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

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

The result of the computation.

           X Matrix

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

             0.9996

             2.0002

             1.0002

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

      Number of iterations : 25

      >>