การคำนวณเชิงตัวเลข (๒๕) ตัวอย่างการแก้ปัญหา สมการพีชคณิตไม่เชิงเส้นด้วย Function fzero ของ GNU Octave MO Memoir : Tuesday 4 February 2563

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

  

รูปที่ ๑ ข้อจำกัดของฟังก์ชัน fzero ของโปรแกรม MATLAB ที่เขียนไว้ในเว็บ MathWorks

งานวิทยานิพนธ์สมัยที่ผมทำนั้นเขียนด้วยภาษา FORTRAN 77 เวลาที่จะแก้สมการหรือระบบสมการก็ยังต้องเขียนอัลกอริทึมกันเอง จะโชคดีหน่อยที่ทางมหาวิทยาลัยซื้อ NAG Library มาให้ใช้ ทำให้การแก้ปัญหาบางอย่างเช่นการแก้ปัญหาระบบสมการเชิงเส้น (system of linear equations) ที่เกี่ยวข้องกับเมทริกซ์หรือการหาคำตอบของสมการไม่เชิงเส้น (nonlinear equation) ก็ไม่จำเป็นต้องเขียนอัลกอริทึมส่วนนั้นเอง แต่ต้องรู้ว่ารูปแบบของสมการที่จะแก้นั้นเป็นอย่างไร จะได้เรียกใช้ libraray ได้ถูกตัว เพราะถ้าเรียกใช้ตัวที่ไม่เหมาะสม ก็อาจมีปัญหาตามมาได้หลายอย่าง เช่นหน่วยความจำไม่พอ เสียเวลาคำนวณนาน มีปัญหาการความคลาดเคลื่อนเนื่องจากการปัดเศษสะสมมาก หรืออาจถึงขั้นหาคำตอบไม่เจอทั้ง ๆ ที่ในความเป็นจริงนั้นคำตอบมันมี ดังนั้นการเรียนการคำนวณเชิงตัวเลขที่เรียนมา จึงเน้นไปที่การทำงาน จุดอ่อน จุดแข็ง และข้อจำกัดของระเบียบวิธีการแก้ปัญหาแต่ละวิธี เพราะถ้ามันมีวิธีที่ดีที่สุดใช้ได้กับทุกกรณีอยู่เพียงวิธีเดียว มันก็คงไม่ต้องเรียนกันเยอะ เรียนกันเพียงแค่วิธีเดียวก็พอ แต่ปัจจุบันนี้ดูเหมือนว่าจะเปลี่ยนไปเยอะ มีความพยายามที่จะให้สอนนิสิตให้รู้จักใช้โปรแกรมสำเร็จรูป โดยที่ไม่จำเป็นต้องรู้ว่าคำสั่งแต่ละคำสั่งของโปรแกรมนั้นมันทำงานอย่างไร มันมีข้อจำกัดในการใช้อย่างไร ซึ่งผมก็ไม่เห็นด้วย ก็เลยขอยกตัวอย่างบางตัวอย่างมาให้ดูกัน
  

รูปที่ ๒ เทคนิคที่ฟังก์ชัน fzero ของ MATLAB ใช้ ซึ่งอันที่จริงมันก็คืออัลกอริทึมที่เขียนไว้ตั้งแต่สมัยยุคทศตวรรษ ๑๙๖๐ ด้วยภาษา ALGOL 60 เพียงแต่มีการนำมาเขียนใหม่ด้วยภาษาใหม่เท่านั้นเอง ก็ไม่ต่างอะไรจากเอาหนังสือที่เขียนไว้ในภาษาหนึ่งมาแปลเป็นอีกภาษาหนึ่ง หนังสือที่แลปมาก็ยังมีเนื้อหาคงเดิม (อาจมีแก้พิมพ์ผิดบ้าง) แต่โดยรวมไม่ต่างไปจากของเดิม

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

การหาคำตอบของสมการที่ไม่เป็นเชิงเส้น (คือมีสมการเดียว) นั้นมีอยู่ด้วยกันหลายวิธี แบบที่เรียกว่าพื้นฐานการทำงานแตกต่างกันเลย ตรงนี้จะแตกต่างไปจากกรณีของการแก้ปัญหาระบบสมการที่ไม่เป็นเชิงเส้น (คือมีตั้งแต่ ๒ สมการขึ้นไป) ที่ดูเหมือนว่าจะมีหลายวิธี แต่เอาเข้าจริง ๆ มันพอจะจัดออกเป็นสองกลุ่ม กลุ่มแรกมันอิงอยู่บนวิธีของ Newton ที่ต้องมีการใช้ Jabobian matrix วิธีการกลุ่มนี้ที่มันแตกย่อยออกไปมักจะอยู่ตรงที่รายละเอียดปลีกย่อยข้างใน ส่วนอีกกลุ่มหนึ่งนั้นก็หลีกไปเล่นทางเทคนิคที่ใช้กับการหาค่าต่ำสุดของฟังก์ชันเลย แต่ไม่ว่าจะเป็นวิธีการแบบใดก็ตามสำหรับระบบกี่สมการก็ตาม มันไม่มีวิธีการใดเลยที่กล้าพูดได้เต็มปากเลยว่าไม่ว่าจะเป็นโจทย์แบบใดก็ตาม ก็สามารถหาคำตอบเจอและหาได้เจอทุกคำตอบ
  

คำสั่งที่ยกมาเป็นตัวอย่างในวันนี้คือคำสั่ง fzero ที่ใช้สำหรับหารากของสมการไม่เชิงเส้นสำหรับกรณีของระบบ ๑ สมการ ๑ ตัวแปร รายละเอียดของคำสั่งนี้ของโปรแกรม MATLAB แสดงไว้ในรูปที่ ๑ และ ๒ แต่โปรแกรมที่จะเล่นให้ดูในวันนี้ไม่ใช้ MATLAB แต่เป็น GNU Octave ที่เป็น free software ที่ใคร ๆ ก็สามารถดาวน์โหลดมาใช้งานได้ฟรี ๆ ชื่อ Octave นี้ก็นำมาจากชื่อ Prof. Octave Levenspiel ที่เป็นศาสตราจารย์ทางด้านวิศวกรรมเคมีที่เขียนตำรา Chemical Reaction Engineering (ผมเองก็เรียนวิชานี้โดยใช้ตำราที่แกเขียน หนังสือเล่มนี้ก็ยังคงเก็บเอาไว้อยู่) และโปรแกรมนี้ตอนแรกก็กะให้เป็นโปรแกรมให้ผู้เรียนวิศวกรรมเคมีใช้กันโดยไม่ต้องไปจ่ายตังค์ชื้อ MATLAB ที่มีราคาแพง (แต่ทำไมหลายสถาบันกลับชอบจ่ายเงิน) เรื่องราวเกี่ยวกับอาจารย์ท่านนี้เคยเขียนไว้ครั้งหนึ่งในเรื่อง "ทำไมไดโนเสาร์จึงสูญพันธุ์" (Memoir ปีที่ ๕ ฉบับที่ ๕๔๐ วันจันทร์ที่ ๒๖ พฤศจิกายน พ.ศ. ๒๕๕๕)
  

รูปที่ ๓ รูปแบบหน้าตาของโปรแกรม GNU Octave ที่รูปแบบการเขียนและการทำงานนั้นใกล้เคียงกับ MATLAB ที่บางคนบอกว่าสามารถนำเอาคำสั่งที่เขียนบน MATLAB มาใช้กับโปรแกรมนี้ได้โดยแทบไม่ต้องแก้ไขอะไร หรือทำกลับกันก็ได้ ที่สำคัญก็คือมันเป็นโปรแกม Open source ที่สามารถดาวน์โหลดมาใช้งานได้ฟรี ไม่ต้องจ่ายเงินซื้อ

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

แต่วิธีการนี้ก็มีปัญหาคือถ้าหากกราฟนั้นสัมผัสกับแกน x แต่ไม่ข้ามแกน x มันจะมีจุดที่ทำให้ค่า f(x) เป็นศูนย์ได้เช่นกัน เช่นกรณีของกราฟ f(x) = x² รูปที่ ๒ นั้นบอกว่าฟังก์ชัน fzero ของ MATLAB หาจุดดังกล่าวด้วยเทคนิค bisection ร่วมกับ secant และ inverse quadratic interpolation ในบรรดาเทคนิคทั้งสามนี้ เทคนิค bisection จำเป็นต้องมีจุดสองจุดที่อยู่คนละฟากของแกน x ก่อนจึงจะเริ่มค้นหาได้ ด้วยเหตุนี้ในรูปที่ ๑ จึงมีคำเตือนเอาไว้ว่าฟังก์ชัน fzero นั้นอาจหาคำตอบของบางฟังก์ชันไม่ได้ แม้ว่าฟังก์ชันนั้นจะมีคำตอบก็ตาม และก็ได้ยกตัวอย่างสมการ f(x) = x²

  
คู่มือของ GNU Octave (รูปที่ ๔) ไม่ได้ให้รายละเอียดว่าฟังก์ชัน fzero ใช้อัลกอริทึมใดในการหาคำตอบ แต่เดาว่าน่าจะเป็นตัวเดียวกับที่ MATLAB ใช้ (เพราะอัลกอริทึมนี้มีใช้มานานแล้วและเปิดเผยให้ทราบทั่วไป) การทดลองเริ่มจากการให้หารากของฟังก์ชัน f(x) = x² ก่อน (สมการที่ 1 ในรูปที่ ๕) โดยให้เริ่มควานหาคำตอบที่จุด x = 2 ก็พบว่าฟังก์ชัน fzero ให้คำตอบค่า x ที่ทำให้ f(x) มีค่าเป็นศูนย์ก็คือที่ x = 0
  

ต่อมาทดลองให้หารากของฟังก์ชัน f(x) = (x-1)² (สมการที่ 2) ค่า x ที่ทำให้ค่า f(x) ของสมการนี้เท่ากับศูนย์คือ x = 1 แต่พอใฟ้ฟังก์ชัน fzero หาคำตอบโดยเริ่มหาจากจุด x = 2 (สมการที่ 2.1) ก็พบว่ามันไม่สามารถหาคำตอบได้เพราะไม่สามารถหาจุดสองจุดที่อยู่คนละฟากของแกน x ได้ (ตรง error ที่รายงานว่า "not a valid initial bracketing") และพอย้ายจุดเริ่มต้นการคำนวณมาเป็นจุด x = -3 มันก็รายงานว่าไม่สามารถหาคำตอบได้ด้วยเหตุผลเดียวกัน แต่ทำไมในกรณีของสมการ f(x) = x² นั้นมันกลับทำได้ก็ไม่รู้เหมือนกัน
  

รูปที่ ๔ รายละเอียดฟังก์ชัน fzero ของ GNU Octave ที่ให้ไว้แต่วิธีการเขียน แต่ไม่ได้บอกว่าใช้เทคนิคอะไร แต่เท่าที่ลองดูพบว่าน่าจะใช้เทคนิคทำนองเดียวกับที่ MATLAB ใช้

  

รูปที่ ๕ การหาคำตอบของฟังก์ชัน f(x) = x², f(x) = (x-1)², f(x) = x² - 0.0001 และ f(x) = (x-1)² - 0.0001 ด้วยฟังก์ชัน fzero ของ GNU Octave

ทีนี้ลองเปลี่ยนเป็นฟังก์ชัน f(x) = x² - 0.0001 (สมการที่ 3) ที่คราวนี้จะมีจุดตัดแกน x อยู่สองจุดคือที่ x = 0.01 และ -0.01 โดยค่า f(x) จะมีค่าติดลบเมื่อ x อยู่ในช่วงจาก -0.01 ถึง 0.01 และนอกช่วงนี้จะมีค่าเป็นบวกทั้งหมด เมื่อให้ฟังก์ชัน fzero เริ่มต้นหาจากจุด x = 1 ก็ได้คำตอบจุดตัดคือ x = 0.01 (สมการที่ 3.1) และเมื่อเปลี่ยนจุดเริ่มต้นการหาเป็นจุด x = -1 (สมการที่ 3.2) ก็ได้คำตอบออกมาเป็น x = -0.01
  

ที่นี้มาทดลองกับสมการ f(x) = (x-1)² - 0.0001 (สมการที่ 4) ที่คราวนี้จะมีจุดตัดแกน x อยู่สองจุดคือที่ x = 0.99 และ 1.01 โดยค่า f(x) จะมีค่าติดลบเมื่อ x อยู่ในช่วงจาก -0.01 ถึง 0.01 และนอกช่วงนี้จะมีค่าเป็นบวกทั้งหมด เมื่อให้ฟังก์ชัน fzero เริ่มต้นหาจากจุด x = 1 ก็ได้คำตอบจุดตัดคือ x = 0.99 (สมการที่ 4.1) แต่เมื่อเปลี่ยนจุดเริ่มต้นการหาเป็นจุด x = -1 (สมการที่ 4.2) โปรแกรมกลับรายงานว่าหาคำตอบไม่ได้เพราะไม่สามารถหาจุดสองจุดที่อยู่คนละฟากของแกน x ได้ และเมื่อเปลี่ยนค่า x เริ่มต้นการหาเป็นจุด x = 1.1 (สมการที่ 4.3) โปรแกรมก็ยังไม่สามารถหาคำตอบได้ด้วยเหตุผลเดียวกัน

  

รูปที่ ๖ การหาคำตอบของฟังก์ชัน f(x) = 1/x, f(x) = 1/(x-1) และ f(x) = 1/(x-1) - 0.1 ด้วยฟังก์ชัน fzero ของ GNU Octave

ที่นี้เป็นเป็นสมการไฮเปอร์โบลาที่ไม่ตัดแกน x แต่เป็นสมการที่จุดอยู่คนละฟากของแกน x โดยเริ่มจาก f(x) = 1/x ก่อน (สมการที่ 5) โดยให้โปรแกรมเริ่มทดลองหาจากจุด x = 2 (สมการที่ 5.1) ก็พบว่าโปรแกรมให้คำตอบว่าจุดตัดแกนของฟังก์ชันดังกล่าวคือจุด x = 1.8398e-6 ซึ่งในความเป็นจริงนั้นค่า f(x) ที่ตำแหน่งนี้อยู่ห่างแกน x ไปลิบลับเลย การที่มันได้คำตอบนี้มาแสดงว่ามันเริ่มต้นด้วยการหาจุดสองจุดที่อยู่คนละฟากของแกน x ก่อน ซึ่งมันก็หาได้ (ก็เพราะมันมี) จากนั้นก็ตามด้วยวิธี bisection ที่ทำไปเรื่อย ๆ จนไม่สามารถแบ่งได้อีก แต่วิธีการนี้มันใช้ได้ก็ต่อเมื่อฟังก์ชันที่เชื่อมระหว่างจุดสองจุดที่อยู่คนละฟากของแกน x นั้นเป็น "ฟังก์ชันต่อเนื่อง" แต่ในกรณีนี้มันเป็นฟังก์ชันไม่ต่อเนื่อง และคงใช้เงื่อนไขหยุดการคำนวณโดยพิจารณาจากค่า ∆x เพียงอย่างเดียว กล่าวคือถ้าพบว่าเมื่อคำนวณไปเรื่อย ๆ แล้วค่า ∆x มีการเปลี่ยนแปลงไม่เกินขนาดที่กำหนดไว้ก็จะหยุดการคำนวณ แต่ถ้าตรวจสอบค่า f(x) ก็จะรู้เลยว่าโดนฟังก์ชันหลอก เพราะถ้าการคำนวณนั้นนำไปสู่จุดที่กราฟตัดแกนx ก็จะเห็นว่าทั้งขนาดการเปลี่ยนแปลงค่าx และค่า f(x) นั้นเข้าสู่ศูนย์ทั้งคู่
  

สมการที่ (6) เป็นกรณีของฟังก์ชัน y = 1/(x-1) ที่พบว่าเมื่อให้หาเริ่มต้นหาคำตอบที่จุด x = 0.1 (สมการที่ 6.1) โปรแกรมให้คำตอบจุดตัดออกมาที่ x = -9.5932e-17 และเมื่อเปลี่ยนจุดเริ่มต้นเป็น x = 3 (สมการที่ 6.2) ก็ได้คำตอบออกมาเป็น x = 1.0000 ซึ่งผิดทั้งคู่ แม้ว่าตัวโปรแกรมจะมีคำเตือนออกมาก็ตาม
  

ที่นี้ลองปรับสมการนิดนึงเป็นกรณีของฟังก์ชัน y = 1/(x-1) - 0.1 (สมการที่ 7) พอให้เริ่มต้นหาคำตอบที่จุด x = 1 (สมการที่ 6.1) หรือ x = 0 (สมการที่ 6.3) ก็พบว่าโปรแกรมให้คำตอบออกมาเป็น x = 1.0312e-16 cและ -4.3521e-16 ตามลำดับ ซึ่งผิดทั้งคู่ แต่พอให้เริ่มต้นหาคำตอบที่จุด x = 12 (สมการที่ 7.2) ก็สามารถหาคำตอบที่ถูกต้องคือ x = 10

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

แต่สิ่งที่คิดว่าเป็นปัญหาสำคัญในตอนนี้คือ ผู้ที่ทำวิจัยโดยใช้โปรแกรมสำเร็จรูปแก้ปัญหาระบบสมการนั้น ไม่มีความรู้เลยว่าควรจะเลือกใช้ฟังก์ชันไหน (เผลอ ๆ ไม่รู้เลยว่ามันมีตัวเลือกอยู่) และสำหรับปัญหาของตัวเองนั้นโปรแกรมดังกล่าวนั้นมีฟังก์ชันที่เหมาะสมหรือไม่ เพราะถ้ามันไม่มีก็ต้องเขียนขึ้นเอง ที่ผ่านมาเห็นอาศัยว่าเทียบกับบทความที่มีตีพิมพ์ว่าคนอื่นเขาใช้ฟังก์ชันไหนแล้วก็ใช้ตาม โดยที่ไม่รู้เลยว่าทำไปต้องใช้ฟังก์ชันนั้น ผมเองตอนสอบวิทยานิพนธ์ของนิสิตรายหนึ่งถึงกับต้องกล่าวว่า

"paper นั้นเขาใช้คำสั่งนั้นไม่ผิดหรอก แต่ของคุณนั้นใช้ไม่ได้ เพราะพฤติกรรมของสมการมันไม่เหมือนกัน สมการของคุณไม่เหมือนกับของเขา มันแค่คล้ายกัน"