Scilab for Science
วันพฤหัสบดีที่ 28 พฤษภาคม พ.ศ. 2563
วันศุกร์ที่ 15 พฤษภาคม พ.ศ. 2563
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
-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)
สร้างเว็คเตอร์ 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).
สมัครสมาชิก:
บทความ (Atom)