การคำนวณเชิงตัวเลข (๑๒) LU decomposition ร่วมกับ Iterative improvement MO Memoir : Saturday 22 August 2558

Memoir ฉบับนี้จัดทำขึ้นเพื่อขยายความในเอกสารประกอบคำสอนวิชาคณิตศาสตร์สำหรับวิศวกรรมเคมี ในส่วนของการคำนวณเชิงตัวเลข เรื่องการแก้ปัญหาระบบสมการพีชคณิตเชิงเส้น

ในการแก้ปัญหาระบบสมการพีชคณิต A.x = b ด้วยเทคนิคการกำจัดของเกาส์ (Gauss elimination) นั้น จำนวนครั้งของการคูณในขั้นนตอนของการกำจัดไปข้างหน้า (Forward elimination) จะแปรผันตามจำนวนสมการยกกำลัง 2 ในขณะที่จำนวนครั้งของการคูณในขั้นตอนการแทนค่าย้อนกลับ (Back substitution) แปรผันตามจำนวนสมการยกกำลัง 1 หรือกล่าวอีกนัยหนึ่งคืองานส่วนใหญ่ไปหนักอยู่ที่ขั้นตอนการกำจัดไปข้างหน้า
  

ในกรณีที่เราต้องแก้ปัญหาระบบสมการพีชคณิต A.x = b หลายครั้ง โดยแต่ละครั้งนั้นมีเมทริกซ์ A ที่เหมือนกัน แต่แตกต่างกันตรงเวกเตอร์ b แม้ว่าเราจะสามารถใช้เทคนิคการกำจัดของเกาส์ในการแก้ปัญหาได้ แต่จะพบว่ามีการทำงานที่ซ้ำซ้อนกันอยู่ คือขั้นตอนการกำจัดไปข้างหน้าที่เราต้องทำการแปลงเมทริกซ์ A เพื่อให้ค่าตัวเลขในซีกซ้ายล่างของเมทริกซ์ A (พวกที่อยู่ใต้แนวเส้นทแยงมุม) มีค่าเป็น 0 และที่สำคัญคือขั้นตอนนี้เป็นขั้นตอนที่ใช้เวลาคำนวณมากที่สุด
 

ด้วยเหตุนี้จึงได้มีการหาวิธีที่จะทำการเก็บค่าตัวเลขที่ใช้ในการแปลงตัวเลขในตำแหน่งต่าง ๆ ที่อยู่ใต้เส้นทแยงมุมของเมทริกซ์ A ให้มีค่าเป็น 0 เพื่อที่จะได้สามารถนำไปใช้งานในครั้งต่อไปได้เลยโดยไม่ต้องทำการคำนวณใหม่ นั่นคือการแตกเมทริกซ์ A ให้กลายเป็นผลคูณของเมทริกซ์ 2 เมทริกซ์คือ L กับ U หรือ A = L.U

 
จะว่าไปแล้วการแปลงเมทริกซ์ A ให้กลายเป็นผลคูณของเมทริกซ์ 2 เมทริกซ์คือ L กับ U นั้นก็คือการทำการกำจัดไปข้างหน้าในระหว่างการใช้เทคนิคการกำจัดของเกาส์ในการแก้ปัญหา เพียงแต่เราเอาตัวเลขที่เราใช้ในการกำจัดให้ตัวเลขที่อยู่ ใต้แนวเส้นทแยงมุมให้มีค่าเป็น 0 ไปเก็บไว้ในเมทริกซ์ L ส่วนเมทริกซ์ U ก็คือเมทริกซ์ A ที่เหลือจากการกำจัดไปข้างหน้านั่นเอง
  

ดังนั้นจากโจทย์ของเรา A.x = b

จะแปลงเป็น L.(U.x) = b

กำหนดให้ U.x = y

ทำการแก้ระบบสมการ L.y = b เพื่อหาค่า y ก่อน

จากนั้นแก้ระบบสมการ U.x = y เพื่อให้ได้ค่า x

เพื่อให้ได้คำตอบที่ออกมาถูกต้อง (เท่าที่ความอุปกรณ์ที่ใช้ในการคำนวณจะให้ได้) การแก้ปัญหาด้วยเทคนิคการกำจัดของเกาส์นั้นจะเริ่มจากการทำการกำจัดไปข้างหน้าพร้อมกับการหาตัวหลัก (pivoting) ไปด้วยพร้อมกัน และในขณะเดียวกันก็จะทำการแตกเมทริกซ์ A ให้กลายเป็นผลคูณของเมทริกซ์ 2 เมทริกซ์คือ L กับ U เมื่อได้คำตอบ (คือ x) แล้วจะนำคำตอบที่ได้ไปแทนค่าในโจทย์เริ่มต้นเพื่อตรวจสอบว่าผลการคำนวณมีความคลาดเคลื่อนสะสมมากน้อยเท่าใด (ดูผลต่างระหว่างค่า b ที่โจทย์กำหนดมาให้กับค่าที่คำนวณได้จากการแทน x เข้าไปในสมการ) ถ้าพบว่าค่า b ที่ได้จากการแทนค่า x ที่คำนวณได้แตกต่างจากค่า b ที่โจทย์กำหนดมาให้ นั่นแสดงว่าค่า x ที่คำนวณได้นั้นมีความคลาดเคลื่อนสะสมอยู่ ก็จะทำการปรับแก้ด้วยการใช้เทคนิค iterative improvement
 

การปรับแก้คำตอบด้วยเทคนิค iterative improvement นั้นเป็นการแก้ปัญหาระบบสมการ A.dx = db เพื่อหาว่าค่าความคลาดเคลื่อนของคำตอบ (dx) นั้นมีค่าเท่าใด เพื่อที่จะได้ไปหักค่าความคลาดเคลื่อนนั้นออกจากค่าที่คำนวณได้ครั้งแรก กระบวนการปรับแก้ดังกล่าวเป็นกระบวนการที่ทำซ้ำไปเรื่อย ๆ จนกว่าจะพบว่าไม่สามารถปรับลดความคลาดเคลื่อนให้ต่ำลงไปได้อีกแล้ว (เกินความสามารถของเครื่องคำนวณ) และกระบวนการนี้ก็เป็นการทำให้การแก้ปัญหาระบบสมการ A.x = b ของเราเพียงข้อเดียว กลายเป็นการแก้ปัญหาระบบสมการ A.x = b หลายครั้ง โดยที่แต่ละครั้งนั้นมีค่าเวกเตอร์ b ที่แตกต่างกัน
  

เพื่อให้เห็นภาพเรามาลองดูตัวอย่างโดยใช้ระบบสมการพีชคณิต 5 สมการ 5 ตัวแปรที่เคยยกมาก่อนหน้านี้

0.001x₁ + 200x₂ + x₃ - 10x₄ - 200x₅ = 5     สมการที่ (1)

0.6x₁ + 5x₂ + 11.5x₃ + 4x₄ + 3x₅ = -3      สมการที่ (2)

-11x₁ + 3x₂ + 0.8x₃ + x₄ - 15.8x₅ = 7      สมการที่ (3)

2x₁ + 2x₂ + 2x₃ + 0.000067x₄ + 7x₅ = 6     สมการที่ (4)

-3x₁ - x₂ + 2x₃ + 8x₄ + 1.3x₅ = -11     สมการที่ (5)

จะทำการหาคำตอบของระบบสมการ (ค่า x₁, x₂, x₃, x₄ และ x₅) เปรียบเทียบกันระหว่างการใช้ Gauss elimination โดยไม่มีการหาตัวหลัก จากนั้นจะใช้เทคนิค Iterative improvement เพื่อทำการปรับแก้คำตอบ การคำนวณกระทำโดยใช้โปรแกรม Spreadsheet ของ OpenOffice 3.3.0
  

รูปที่ ๑ เป็นผลการคำนวณโดยไม่มีการหาตัวหลัก เมื่อนำค่า x₁, x₂, x₃, x₄ และ x₅ ที่คำนวณได้แทนค่ากลับเข้าไปในสมการที่ (1)-(5) พบว่าค่า b นั้นมีความคลาดเคลื่อนอยู่ (ตัวเลขสีแดงที่มุมขวาล่างของรูป) ดังนั้นจึงทำการแก้ระบบสมการ L.U.dx = db เพื่อหาค่า dx (ณ ขณะนี้เราไม่จำเป็นต้องทำการแตกเมทริกซ์ A ให้กลายเป็นผลคูณของเมทริกซ์ L กับ U เพราะเราได้ทำไปแล้วในระหว่างการหาคำตอบ x ชุดแรก) รูปที่ ๒ แสดงขั้นตอนคำนวณเพื่อคำนวณหาค่า dx ครั้งแรกเพื่อนำไปปรับคำตอบ x ที่ได้จากการคำนวณในครั้งแรก และเมื่อได้ค่า x ชุดใหม่แล้วก็ทำการทดสอบว่ายังมีความคลาดเคลื่อนอยู่หรือไม่ คือดูว่าค่า db เป็นศูนย์หรือไม่ ซึ่งก็พบว่ายังมีความคลาดเคลื่อนอยู่ (ค่า db ยังไม่เป็นศูนย์)
  

รูปที่ ๓ แสดงการทำ Iterative improvement ซ้ำเป็นครั้งที่สอง เพื่อหาว่าค่า x ที่ได้มาจากการปรับแก้ครั้งแรกนั้นยังมีความคลาดเคลื่อนหลงเหลืออยู่เท่าใด ซึ่งก็พบว่ามีค่าความคลาดเคลื่อนลดลง (จาก 10^-11-10^-12 ลงเหลือ 10^-16-10^-17) และเมื่อนำค่าคลาดคลาดเคลื่อนที่คำนวณได้ใหม่นี้ไปหักออกจากค่า x ที่ได้จากการปรับแก้ครั้งแรก ก็จะได้ค่า x ที่ผ่านการปรับแก้สองครั้ง แต่เมื่อนำค่า x ที่ได้จากการปรับแก้ครั้งที่สองนี้แทนค่ากลับเข้าไปใหม่กลับพบว่าค่าความคลาดเคลื่อน db นั้นไม่ลดต่ำลง ทั้งนี้เป็นเพราะค่า dx ที่ได้จากการคำนวณครั้งที่สองมีขนาดประมาณ 10^-16 เท่าของค่า x ที่ได้จากการปรับแก้ครั้งแรก และการคำนวณนั้นกระทำที่นัยสำคัญ 16 ตำแหน่ง ดังนั้นเมื่อนำค่า dx ที่ได้จากการคำนวณครั้งที่สองไปหักออกจากค่า x ที่ได้จาก การคำนวณครั้งแรก ผลลัพธ์ที่ได้จึงเป็นค่า x ชุดเดิม กล่าวอีกนัยหนึ่งคือเราไม่สามารถปรับคำตอบให้ถูกต้องได้มากกว่านี้อีกแล้วเนื่องด้วยข้อจำกัดของ machine precision ดังนั้นการคำนวณจะยุติลงเพียงแค่นี้

ไฟล์ที่แนบมาด้วยเป็นไฟล์ Spreadsheet ที่ใช้ในการคำนวณ
 

การคำนวณเชิงตัวเลข (๑๑) เปรียบเทียบ Gauss elimination ที่มีและไม่มีการทำ Pivoting MO Memoir : Friday 21 August 2558

Memoir ฉบับนี้จัดทำขึ้นเพื่อขยายความในเอกสารประกอบคำสอนวิชาคณิตศาสตร์สำหรับวิศวกรรมเคมี ในส่วนของการคำนวณเชิงตัวเลข เรื่องการแก้ปัญหาระบบสมการพีชคณิตเชิงเส้น

เทคนิคการกำจัดของเกาส์ (หรือ Gauss elimination) เป็นเทคนิคหลักในการหาคำตอบของระบบสมการพีชคณิตที่มีจำนวนสมการไม่มาก (สมการอยู่ในรูปแบบ A.x = b) เทคนิคนี้ประกอบด้วยขั้นตอนหลัก ๒ ขั้นตอนคือ ขั้นตอนแรกเป็นขั้นตอนของการกำจัดไปข้างหน้า (Forward elimination) ซึ่งเป็นการทำให้ตัวเลขในซีกซ้ายล่างของเมทริกซ์ A มีค่าเป็น 0 และขั้นตอนที่สองคือการแทนค่าย้อนกลับ (Back substitution) ที่เป็นการหาค่า ย้อนกลับจากค่า x_n มาจนถึง x₁ 
   

แต่การแก้โจทย์นั้นไม่ใช่ว่าโจทย์เรียงลำดับสมการมาอย่างไรเราก็แก้โจทย์ไปตามนั้น แต่ควรมีการตรวจสอบและจัดเรียงลำดับสมการใหม่ (ถ้าจำเป็น) ก่อนที่จะทำกระบวนการการกำจัดไปข้างหน้าในแต่ละครั้ง ที่เรียกว่าการหาตัวหลัก (Pivoting) ทั้งนี้เพื่อให้ตัวเลขที่อยู่ในแนวทแยงมุมของเมทริกซ์ A มีขนาดใหญ่ที่สุด (ไม่คิดเครื่องหมาย) ทั้งนี้เพราะตัวเลขที่อยู่ในตำแหน่งทแยงมุมของเมทริกซ์ A จะทำหน้าที่เป็นตัวหารในระหว่างการทำการกำจัดไปข้างหน้า การที่ตัวหารมีขนาดที่ใหญ่กว่าตัวตั้งจะช่วยลดความคลาดเคลื่อนในการคำนวณให้น้อยลง
  

เพื่อให้เห็นภาพเรามาลองดูตัวอย่างกันดีกว่า โดยสมมุติว่าเรามีระบบสมการพีชคณิต 5 สมการ 5 ตัวแปรดังต่อไปนี้

0.001x₁ + 200x₂ + x₃ - 10x₄ - 200x₅ = 5      สมการที่ (1)

0.6x₁ + 5x₂ + 11.5x₃ + 4x₄ + 3x₅ = -3       สมการที่ (2)

-11x₁ + 3x₂ + 0.8x₃ + x₄ - 15.8x₅ = 7      สมการที่ (3)

2x₁ + 2x₂ + 2x₃ + 0.000067x₄ + 7x₅ = 6      สมการที่ (4)

-3x₁ - x₂ + 2x₃ + 8x₄ + 1.3x₅ = -11      สมการที่ (5)

จะทำการหาคำตอบของระบบสมการ (ค่า x₁, x₂, x₃, x₄ และ x₅) เปรียบเทียบกันระหว่างการใช้ Gauss elimination โดยไม่มีการหาตัวหลัก และการใช้ Gauss elimination ร่วมกับการทำการหาตัวหลักแบบบางส่วน (Partial pivoting) การคำนวณกระทำโดยใช้โปรแกรม Spreadsheet ของ OpenOffice 3.3.0

รูปที่ ๑ เป็นการหาคำตอบด้วยเทคนิค Gauss elimination โดยไม่มีการหาตัวหลัก ตัวเลขสีเขียวที่อยู่ข้างหน้าสุดคือตัวเลขที่ใช้ในการกำจัดตัวแปร x_i (i = 2-5) ออกทีละตัวจากแต่ละสมการ ตัวเลขสีแดงคือตัวเลขที่อยู่ในแนวทแยงมุมที่ทำหน้าที่เป็นตัวหาร จะเห็นว่าเมื่อนำค่า x_i (i = 1-5) ที่คำนวณได้แทนค่ากลับเข้าไปในโจทย์ พบว่าค่าที่ได้จากสมการที่ (2)-(5) นั้นมีความคลาดเคลื่อนไปจากค่าที่ควรจะเป็น (ค่า df ซึ่งเป็นผลต่างระหว่างค่าที่ได้จากการนำเอาค่า x₁, x₂, x₃, x₄ และ x₅ ที่คำนวณได้ แทนกลับเข้าไปในโจทย์ กับค่าที่โจทย์กำหนดมาให้คือค่าทางด้านขวามือของเครื่องหมายเท่ากับ) ทั้งนี้เป็นผลของค่าความคลาดเคลื่อนที่สะสมอยู่ในตัวแปร x แต่ละตัว
  

รูปที่ ๒ เป็นการหาคำตอบด้วยเทคนิค Gauss elimination ร่วมกับการหาตัวหลักแบบบางส่วน ในที่นี้พบว่าก่อนที่จะเริ่มทำการกำจัดไปข้างหน้าเป็นครั้งแรก ควรที่จะทำการจัดลำดับสมการใหม่ก่อน โดยเลื่อนสมการที่ (3) ให้มาเป็นสมการที่ (1) ทั้งนี้เพราะค่าสัมประสิทธิ์หน้าตัวแปร x₁ ในสมการที่ (3) ที่มีขนาดเท่ากับ 11 นั้นมีขนาดที่ใหญ่กว่าขนาดของค่าสัมประสิทธิ์หน้าตัวแปร x₁ ในสมการที่ (1) ที่มาค่าเพียง 0.001
  

และหลังจากกำจัด x₁ ออกจากสมการที่ (2) ถึง (5) ไปแล้ว (ตัวเลขในหลักที่ 1 ของแถวที่ 2-5 ของเมทริกซ์ A มีค่าเป็น 0) ก่อนที่จะทำการกำจัดไปข้างหน้าในครั้งต่อไป (เพื่อกำจัดตัวแปร x₂ ออกจากสมการที่ (3) ถึง (5)) ก็ต้องทำการตรวจสอบว่าค่าสัมประสิทธิ์ของตัวแปร x₂ ในสมการที่เหลือ ของสมการใดมีขนาดที่ใหญ่ที่สุด โดยในที่นี้พบว่าค่าสัมประสิทธิ์ของตัวแปร x₂ ของสมการที่ (2) มีขนาดที่ใหญ่ที่สุดเมื่อเทียบกับของสมการอื่นที่เหลือ ดังนั้นจึงไม่จำเป็นต้องมีการจัดลำดับสมการใหม่ สามารถทำการกำจัดไปข้างหน้าเป็นครั้งที่สองได้เลยเพื่อกำจัดตัวแปร x₂ ออกจากสมการที่ (3) ถึง (5)
  

ทำนองเดียวกันก่อนที่จะกำจัดตัวแปร x₃ ออกจากสมการที่ (4) ถึง (5) ก็ต้องทำการตรวจสอบว่าค่าสัมประสิทธิ์ของตัวแปร x₃ ในสมการที่เหลือ ของสมการใดมีขนาดที่ใหญ่ที่สุด โดยในที่นี้พบว่าค่าสัมประสิทธิ์ของตัวแปร x₃ ของสมการที่ (3) มีขนาดที่ใหญ่ที่สุดเมื่อเทียบกับของสมการอื่นที่เหลือ ดังนั้นจึงไม่จำเป็นต้องมีการจัดลำดับสมการใหม่ สามารถทำการกำจัดไปข้างหน้าเป็นครั้งที่สองได้เลยเพื่อกำจัดตัวแปร x₃ ออกจากสมการที่ (4) และ (5)
  

ในขั้นตอนสุดท้ายซึ่งเป็นการกำจัดตัวแปร x₄ ออกจากสมการที่ (5) ในที่นี้พบว่าสัมประสิทธิ์หน้าตัวแปร x₄ ของสมการที่ (5) (ตัวเลขสีน้ำเงิน) มีขนาดใหญ่กว่าของสัมประสิทธิ์หน้าตัวแปร x₄ ในสมการที่ (4) ดังนั้นก่อนจะทำการกำจัดไปข้างหน้าจึงทำการจัดลำดับสมการใหม่ โดยสลับลำดับสมการที่ (4) และ (5)
 

รูปที่ ๒ (ต่อ) เป็นการแทนค่าย้อนกลับเพื่อหาค่าตัวแปร x₁, x₂, x₃, x₄ และ x₅ และนำคำตอบที่ได้แทนค่ากลับเข้าไปในโจทย์ ในที่นี้พบว่าค่าที่ได้จากทุกสมการนั้นตรงกับค่าที่โจทย์กำหนดให้ พึงสังเกตว่าในกรณีของตัวอย่างที่ยกมานี้ความแตกต่างของ x₁, x₂, x₃, x₄ และ x₅ ที่ได้จากการคำนวณโดยมีและไม่มีการหาตัวหลักนั้นอยู่ที่ทศนิยม ๓ ตำแหน่งสุดท้าย

ไฟล์ที่แนบมาด้วยเป็นไฟล์ Spreadsheet ที่ใช้ในการคำนวณ
  

รูปที่ ๑ การหาคำตอบด้วยเทคนิค Gauss elimination โดยไม่มีการหาตัวหลัก
  

รูปที่ ๒ การหาคำตอบด้วยเทคนิค Gauss elimination ร่วมกับการหาตัวหลัก ในที่นี้มีการสลับลำดับสมการสองครั้ง คือการสลับลำดับครั้งแรกกระทำก่อนที่จะทำการกำจัดไปข้างหน้าครั้งแรก โดยได้ทำการย้ายสมการที่ (3) ขึ้นไปแทนที่สมการที่ (1) และดันสมการที่ (1) ลงมาเป็นสมการที่ (2) และดันสมการที่ (2) ลงมาเป็นสมการที่ (3) การสลับลำดับสมการครั้งที่สองเกิดขึ้นอีกครั้งก่อนที่จะทำการกำจัดไปข้างหน้าในครั้งที่สี่ โดยทำการสลับลำดับสมการที่ (4) และ (5) ทั้งนี้เพราะค่าสัมประสิทธิ์ของตัวแปร x₄ ในสมการที่ (5) (ตัวเลขสีน้ำเงิน) มีขนาดใหญ่กว่าค่าสัมประสิทธิ์ของตัวแปร x₄ ในสมการที่ (4)
  

รูปที่ ๒ (ต่อ) การหาคำตอบด้วยเทคนิค Gauss elimination ร่วมกับการหาตัวหลัก

ดาวน์โหลดไฟล์ pdf ที่ link นี้