Temp

Scilab for Computational Science

เนื้อหาเกี่ยวข้องกับ ไซแลบ  ตั้งแต่การติดตั้งซอพท์แวร์ วิธีการใช้งาน รวมทั้งตัวอย่างการประยุกต์ใช้ทางวิทยาศาสตร์ในสาขาต่างๆ

การติวไซแลบเพื่อวิทยาศาสตร์การคำนวณ

 
ถอดความเอกสารนี้โดย หัสชัย สิทธิรักษ์ บล็อกถัดไปจะแสดงเมนูของเอกสารนี้

บทที่ 2 การคำนวณเททริกซ์

บทที่ 2

การคำนวณเมทริกซ์

2.1 บทนำ

            ในตอนนี้ออกแบบมาเพื่อให้ผู้ใช้ได้คุ้นเคยกับการคำนวณเมทริกซ์ในไซแลบ.

2.2 เมทริกซ์ และเว็คเตอร์

            แม้ว่าไซแลบจะเป็นเครื่องคิดเลขที่มีประโยชน์แล้วก็ตาม แต่ประสิทธิ์ภาพหลักคือให้แนวทางที่ง่ายในการทำงานกับเมทริกซ์และเว็คเตอร์  .ในบทที่แล้วได้เห็นการใช้เว็คเตอร์ในการพล็อตกราฟแล้ว

 จำการพิมพ์คำสั่งตามที่ประกฏให้เห็นที่นี่ สังเกตและเข้าใจการตอบสนองของไซแลบ  ถ้าหากพิมพ์พลาดผิดไป ก็ง่ายในการแก้ไขโดยใช้คีย์ลูกศร และคีย์ [Del]  และให้ทดลองใช้คำสั่งขอความช่วยเหลือเพื่อหาว่าคำสั่งที่ไม่คุ้นเคยเป็นอย่างไร นอกจากนี้อาจจะถามจากผู้รู้

            2.2.1 การแก้สมการ Solving Equations.

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

 x1 + 2x2 - x3 = 1

-2x1 - 6x2 + 4x3 = -2

-x1 - 3x2 + 3x3 = 1

            แม้ว่าปัญหาเหล่านี้สามารถที่จะแก้ได้ด้วยเมือ โดยการขจัดตัวไม่ทราบออกไป  นั้นน่าเบื่อแล้วข้อผิดพลาดมักเกิดขึ้นเสมอ.

            ในการเรียนคณิตศาสตร์ในปีแรกปัญหาที่เขียนขึ้นตามข้อกำหนด เมทริกซ์-เว็คเตอร์ โดยนำเข้าสู่เมทริกซ์ a  และเว็คเตอร์ bดังนี้

    a = 1    2    -1                  b =   1
         -2   -6    4                          -2

         -1   -3    3                          1

 

            จากนี้ต้องการหาคำตอบเว็คเตอร์ x = [x1, x2, x3]  ดังนั้นจะได้ตามสมการ

            Ax = b.

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

// Set up a system

A = [ 1 2 -1; -2 -6  4 ; -1 -3  3 ] // of equations.

b = [ 1; -2; 1 ]

x = A\b                                                             // Find x with A x = b.

            ไซแลบควรจะให้คำตอบคือ

            ความคิดที่สอดคล้องกันจากการคำนวณด้วยมือคือการแทนค่าคำตอบจากการคำนวณเข้าไปในสมาการ และดูว่าทุกสมการถูกต้อง  ในไซแลบสามารถทำได้โดยตรวจสอบผลคูณเว็คเตอร์เมทริกซ์ Ax ที่จะให้ผลลัพธ์เป็นเวคเตอร์ b,  ที่ดีไปกกว่านั้น อาจตรวจสอบดูว่าt b − Ax เท่ากับศูนย์จริงหรือไม่.

A*x                                                                  // ตรวจสอบค่า Ax และ b เท่ากัน

b

 b - Ax                                                               // ตรวจสอบ b - A x = 0 หรือไม่ .

เว็คเตอร์ b − Ax รู้จักกันว่าเป็นส่วนเหลือหรือส่วนเกิน (residue)

            2.2.2 เมทริกซ์และเว็คเตอร์.

            การคำนวณทางวิทยาศาสตร์จำนวนมากเกี่ยวข้องกับการแก้ปัญหาระบบสมการขนาดใหญ่ ที่มีตัวไม่ทราบค่านับล้านๆ  ในตอนนี้ได้ให้ข้อปฏิบัติกับคำสั่งที่จำเป็นในการทำงานกับเมทริกซ์ที่ใหญ่กว่า  ตัวอย่างชุดของรหัสต่อไปนี้จะกำหนดเมทริกซ์ A ขนาด 4 × 4 และจากนั้นหาคุณสมบัติที่สำคัญบางอย่าง เมื่อ[ ]  นำมาใช้ในการกำหนดเมทริกซ์  ทั้งช่วงว่างหรือคอมม่า (,) สามารถใช้ในการแยกองค์ประกอบแถวที่นำเข้ามาในเมทริกซ์  เครื่องหมายเซมิโคลอน  (;) ใช้สำหรับเริ่มต้นแถวใหม่

A = [1 2 3 4 ; 1 4 9 16 ; 1 8 27 64 ; 1 16 81 256 ]

A’

det(A)

spec(A)

[D, X] = bdiag(A)

inv(A)

            แม้ว่าในกรณีเมทริกขนาด 4x4 ที่ง่าย แต่ก็เป็นงานน่าเบื่อในการหาค่า ดีเทอร์มีแนนซ์ อินเวิร์ส  แต่ในไซแลบสามารถหาค่าได้ง่ายและรวดเร็ว และใช้ดัชนีเพื่อแสดงส่วนต่างๆของเมทริกซ์ขนาดใหญ่ได้   ตัวอย่างเช่น จากการทดลอง

A(2,3), A(1:2,2:4), A(:,2), A(3,:), A(2:$,$)

            โดยทั่วไป A( i:j, k:l )  หมายถึงตำแหน่งระหว่างแถว i ถึงแถว j  และ คอลัมน์ k ถึงคอลัมน์ l  ช่วง (ranges) สามารถแทนที่ได้โดยใช้ ‘:’ ถ้าแถวหรือคอลัมน์ทั้งหมดรวมอยู่ภายใน  สำหรับแถวหรือคอลัมน์สุดท้ายกำหนดด้วยเครื่องหมาย $.

            แบบฝึกหัดที่ 4.

1. ให้แสดงผลมุมล่างซ้ายขนาด  2 × 3 ของ A ได้อย่างไร? หาค่าดีเทอร์มีแนนซ์ชองมุมบนซ้ายบล็อกขนาด 3 × 3 ของ A ได้อย่างไร?

2. จงหาคำสั่งไซแลบเพื่อคำนวณแถวที่เลือกจากเมทริกซ์ (rowchelon) รูปแบบแถวที่เลือก (echelon) ไม่ได้มีประโยชน์ทางปฏิบัติมาก่อน แต่ก็ง่ายที่จะตรวจสอบการบ้านในวิชาอื่นๆ

            2.2.3 การสร้างเมทริกซ์

            ไซแลบได้ให้แนวทางที่รวดเร็วในการเมทริกซ์พิเศษและเว็คเตอร์ไว้ ดังเช่น

c = ones(4,3)

d = zeros(20,1)

I = eye(5,5)

D = diag( [2 1 0 -1 -2] )

L = diag( [1 2 3 4], -1 )

U = rand(5,5) // U is a 5 x 5 matrix                       // of uniformly distributed random numbers.

R = rand(5,5,’normal’);                                       // R is a 5 x 5 matrix

                                                                        // of normally distributed random numbers.

(ข้อสนเทศเพิ่มเติมสำหรับคำสั่งเหล่านี้สามารถพบได้จากคำสั่ง help.)

            การคูณ เมทริกซ์-เมทริกซ์ และ เมทริกซ์-เว็คเตอร์ เป็นไปตามที่คาดหวัง ให้ผลเป็นมิติที่สอดคล้องกัน  เมทริกซ์และเว็คเตอร์ทั้งสองสามารถคูณกับค่าสเกลล่า  ในตัวอย่างต่อไป B เป็นอีกเมทริกซ์ขนาด 4×4 และ c เป็นเว็คเตอร์คอลัมน์ความยาว 4 (นั่นคือเมทริกซ์หนึ่งขนาด 4x1)

B = [1 1 0 0

0 2 1 0

0 0 3 1

0 0 0 4 ]                                   // Rows can be separated by

// making a new line as well as using ‘;‘.

c = [1; 0; 0; -1]

5*B                                          // Multiply scalars and matrices.

B*c                                           // Multiply matrices and vectors.

A*B                                          // Matrix by matrix multiplication.

B*A                                          // Note: AB is not the same as BA !

            บางคำสั่งต่อไปนี้ใช้ไม่ได้ ไม่อันตรายอะไรที่จะทดลองพิมพ์ดู

c*c // Wrong dimensions for matrix multiplication.

c*A // Wrong dimensions for matrix multiplication.

c’ // The transpose of c, i.e. a row vector.

c’*A // Now matrix multiplication is permitted.

c*c’ // This multiplies a 4 x 1 by a 1 x 4 matrix.

c’*c // What is this quantity usually called?

            แบบฝึกหัดที่ 5.

1. จงคำนวณ inv(A)*A  และ A*inv(A).  ให้ผลลัพธ์ตามที่คาดหวังหรือไม่? (จงใช้คำสั่ง help ถ้าหากไม่สามารถเดาทายได้ว่าคำสั่ง inv ทำอะไร)

2. จงตรวจสอบเมทริกซ์  A และ B ข้างบน ที่ (AB)−1 = B−1A−1.

3. จงตรวจสอบคำสั่งไซแลบ A^(-1) สมารถนำมรใช้เพื่อทำอินเวิรสของ A.

            2.2.4 ระบบของสมการ

            พิจารณาอีกครั้งสำหรับปัญหาการแก้ระบบสมการ Ax = b  อันเป็นแนวทางที่ปรากฏชัดเจนทางคณิตศาสตร์เพื่อคำนวณ A−1b (b/A) . ตามกรณี 3x3 ที่แล้ว นั้นควรจะได้ตรวจสอบคำตอบว่าถูกต้องหรือไม่ ซึ่งทำโดยการคำนวณหาส่วนเหลือ (residual) Ax − b, สำหรับการคำนวณหาคำตอบ x  เนื่องข้อผิดพลาดจากการปัดเลข  ส่วนเหลือนี้อาจจะไม่เท่ากับ 0 ทีเดียว แต่ทุกองค์ประกอบควรจะมีค่าน้อยเท่าที่เป็นได้  กล่าวคือน้อยกว่า 1 × 10−15.

            x = A^(-1)*b                                          // Solve the equation A x = b (A*x = b)

// ตรวจสอบ x ถูกต้อง
//  แน่ใจว่าสอดคล้องตามสมการ.

            resid = b - A*x                                        // ตรวจสอบว่าส่วนเหลือเป็นศูนย์

                                                                        // อย่างน้อยที่สุดภายใต้การปัดเลข.

            ที่จริงแล้วค่อนข้างเร็ว และละเอียดถูกต้องมากขึ้นในการแก้สมการโดยใช้ตัวกระทำ  backslash or ( \ ) ที่นำเสนอไว้ในตอนที่ 1 มากกว่าที่จะคำนวณจากการอินเวิร์ส A−1.  ในบรรทัดต่อไปใช้คำสั่งbackslash ในการแก้สมการ Ax = b,  และตรวจสอบว่าให้คำตอบเดียวกันหรือไม่ดังที่ได้คำนวณจาก A−1.

x1 = A \ b                                   // เป็นแนวทางที่เร็วกว่าในการแก้สมการ A x = b

x - x1                                         // จากนี้คำตอบทั้งสองบเหมือนกัน(จากข้อผิดพลาดการปัดเลข).

            แบบฝึกหัดที่ 6.

            1. ให้แก้ระบบสมการ ดังนี้

            ให้ใช้ตัวกระทำ \.  ให้ตรวจสอบการคำนวณผลเฉลยว่าสอดคล้องตามสมการจริง

(นั่นคือตรวจสอบว่า b − Ax = 0, ต่างไปจากข้อผิดพลาดจากการปัดเลข).

            2. ให้ใช้ rand(n,m,’normal’) เพื่อสร้างเมทริกซ์ขนาด 700 x 700  A,  และเว็คเตอร์คอลัมน์  b ขนาดความยาว 700.  ให้แก้ระบบของสมการ Ax = b โดยใช้คำสั่ง x = A\b  และ x = inv(A)*b, โดยจับเวลาแต่ละคำสั่งโดยใช้คำสั่ง timer()  และคำสั่งใดทำงานได้เร็วกว่าเพราะเหตุใด?

แนะ: ก่อนที่จะทำงานกับเมทริกซ์ขนาดใหญ่ให้ใช้คำสั่ง clear เพื่อเอาตัวแปรทั้งหมดออกจากหน่วยความจำของไซแลบ  เป็นการล้างให้มีหน่วยความจำใช้งานเพียงพอเพื่อไปใช้กับปัญหาหนักๆ ถ้ามีปัญหาอีกก็ให้พยายามลดระบบของสมการให้เล็กลงให้เหมาะกับคอมพิวเตอร์ที่ใช้.

            3. กำหนดให้เมทริกซ์ฮิลเบิร์ทขนาด 10×10 .  องค์ประกอบ ij กำหนดขึ้นโดย 1/(i + j−1). สามารถใช้ฟังก์ชัน feval ทำให้เกิดเมทริกซ์ที่กำหนดฟังก์ชัน

 function z=f(x,y)

     z=1.0/(x+y-1)

 endfunction

สร้างเว็คเตอร์ b ทางด้านขวามือ ซึ่งเป็นคอลัมน์แรกของ A,  b =  A(:,1)  การจะแก้สมการ Ax = b. นี่คือการทดสอบปัญหาทั่วไปในทางพีชคณิตเชิงเส้น (อะไรเป็นผลเฉลยที่แท้จริงของ x ของสมการ Ax = b)

            (a) แก้สมการ Ax = b โดยใช้คำสั่ง backslash และจากนั้นตรวจสอบการคำนวณผลเฉลยว่าสอดคล้องตามสมการโดยการคำนวณส่วนเกิน ( b – Ax)

            (b) จากนี้แก้สมการโดยใช้คำสั่งA^(-1)* ไม่ต้องคำนวณส่วนเกินสำหรับผลเฉลยที่สองนี้.

            (c) อันไหนที่ให้ส่วนเกินที่ดีหว่า (นั่นคือน้อยกว่า)

            (d) วิธีใดให้ค่าผิดพลาดน้อยกว่า

คำอธิบาย: ขณะที่  b เป็นคอลัมน์แรกของ A,  โดยที่สามารถที่จะทำให้ค่าที่แท้จริงของ x โดยปราศจากการคำนวณ  (อธิบายอย่างสังเขป)  ค่าผิดพลาดคือความแตกต่างระหว่างผลเฉลยที่คำนวณของไซแลบและผลเฉลยที่แท้จริง (ผลเฉลยที่แท้ คือเรื่องหนึ่งที่หาออกมาจากความคิดเชิงนามธรรม สามารถที่จะพิมพ์เข้าไปโดยตรงในไซแลบ )  ขนาดของค่าผิดพลาดที่เป็นค่าสูงสุดขององค์ประกอบเว็คเตอร์ค่าผิดพลาด (ด้วยเช่นไซแลบได้ให้ค่าอินเวิร์สเมทริกซ์ของฮิลเบิร์ทด้วยคำสั่ง testmatrix(’hilb’,n).