การคำนวณเชิงตัวเลข (๒๕) ตัวอย่างการแก้ปัญหา สมการพีชคณิตไม่เชิงเส้นด้วย 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 นั้นเขาใช้คำสั่งนั้นไม่ผิดหรอก แต่ของคุณนั้นใช้ไม่ได้ เพราะพฤติกรรมของสมการมันไม่เหมือนกัน สมการของคุณไม่เหมือนกับของเขา มันแค่คล้ายกัน"

การคำนวณเชิงตัวเลข (๒๔) ตัวอย่างการแก้ปัญหา สมการพีชคณิตไม่เชิงเส้นด้วยระเบียบวิธี Müller และ Inverse quadratic interpolation MO Memoir : Sunday 2 February 2563

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

การหาคำตอบของสมการพีชคณิตไม่เชิงเส้น (nonlinear algebraic equation) มีด้วยกันหลากหลายวิธี วิธีการ Secant นั้นก็เป็นวิธีการหนึ่งที่มีการนำไปใช้ในซอร์ฟแวร์สำเร็จรูปที่นิยมให้กันอยู่ในปัจจุบัน วิธีการ Secant นี้เริ่มด้วยการเดาจุดตั้งต้นขึ้นมาก่อน ๒ จุด (จุดที่ ๑ และ ๒) จากนั้นก็คำนวณว่าเส้นตรงที่ลากผ่านสองจุดนี้ไปตัดแกน x ที่ตำแหน่งใด (จุดที่ ๓) จากนั้นคำนวณค่า y ที่จุดนั้น ถ้าพบว่าจุดที่ ๓ ที่ได้มานี้ไม่ใช่คำตอบ ก็โยนจุดที่ ๑ ทิ้งไป แล้วเริ่มต้นคำนวณใหม่ว่าเส้นตรงที่จากผ่านจุดที่ ๒ และ ๓ นั้นไปตัดแกน x ที่ตำแหน่งใด (จุดที่ ๔) ทำอย่างนี้ซ้ำไปเรื่อย ๆ จนกว่าจะพบคำตอบ
    

Müller นำเสนอวิธีใหม่ในปีค.ศ. ๑๙๕๖ (พ.ศ. ๒๔๙๙) วิธีการของ Müller นั้นคล้ายกับวิธีการ Secant เพียงแต่วิธีการ Müller ใช้จุดทั้งสิ้น ๓ จุด จากนั้นสร้างสมการฟังก์ชันพหุนามกำลัง 2 (Quadratic eqaution หรือ Equation of degree 2) ที่ลากผ่านจุด ๓ จุดนั้น แล้วคำนวณว่ารากของสมการกำลัง 2 นั้นอยู่ที่ตำแหน่งใด (คือจุดที่กราฟของสมการกำลัง 2 ตัดแกน x) จากนั้นคำนวณค่า y ที่จุดนั้น ถ้าพบว่าจุดที่ ๔ ที่ได้มานี้ไม่ใช่คำตอบ ก็โยนจุดที่ ๑ ทิ้งไป แล้วเริ่มต้นคำนวณใหม่ว่าเส้นโค้งสมการกำลัง 2 ที่จากผ่านจุดที่ ๒, ๓ และ ๔ นั้นไปตัดแกน x ที่ตำแหน่งใด (จุดที่ ๕) ทำอย่างนี้ซ้ำไปเรื่อย ๆ จนกว่าจะพบคำตอบ รูปที่ ๑ ข้างล่างเป็นการเปรียบเทียบการทำงานระหว่างวิธี Secant และ Müller

  

รูปที่ ๑ เปรียบเทียบการทำงานระหว่างวิธี Secant และ Müller

หลังจากที่เดาจุดเริ่มต้นการคำนวณขึ้นมา ๓ จุด ( x₀, x₁ และ x₂ ) การคำนวณจุดที่สมการกำลังสองไปตัดแกน x (จุด x₃) ตามระเบียบวิธี Müller ทำได้ดังนี้

เมื่อ sign(b) คือเครื่องหมายของตัวแปร b ค่านี้จะเท่ากับ -1 ถ้า b น้อยกว่าศูนย์ และเท่ากับ 1 ถ้า b มากกว่าศูนย์ (อันที่จริงรากของสมการกำลัง 2 นั้นมีสองค่า ขึ้นอยู่กับเครื่องหมายหน้า square root ของตัวหารว่าเป็นบวกหรือลบ แต่ในที่นี้การที่ให้มันมีเครื่องหมายเดียวกับค่า b ก็เพื่อให้ส่วนที่เป็นตัวหารมีขนาดใหญ่ จะได้ลดปัญหาการปัดเศษ)
  

วิธีการหนึ่งที่เป็นญาติกับวิธีการ Müller คือ Inverse quadratic interpolation ความแตกต่างอยู่ตรงที่แทนที่จะใช้จุดข้อมูล ๓ จุดสร้างฟังก์ชันสำหรับคำนวณค่า y ที่ตำแหน่ง x ใด ๆ ก็เปลี่ยนเป็นสร้างฟังก์ชันสำหรับคำนวณค่า x ที่ตำแหน่ง y ใด ๆ แทน โดยตำแหน่งที่กราฟตัดแกน x ก็คือตำแหน่งที่ y เท่ากับศูนย์
  

หลังจากที่เดาจุดเริ่มต้นการคำนวณขึ้นมา ๓ จุด (x₀, f(x₀)), (x₁, f(x₁)) และ (x₂, f(x₂)) การคำนวณจุดที่กราฟตัดแกน x หรือจุด x₃ คำนวณได้จากสมการต่อไปนี้

ถ้าเปรียบเทียบกับรูปแบบสมการกับรูปแบบของ Müller แล้วจะเห็นว่ารูปแบบ Inverse quadratic interpolation มีความซับซ้อนที่น้อยกว่า
  

จุดเด่นข้อหนึ่งของวิธีการของ Müller หรือ Inverse quadratic interpolation คืออัตราการลู่เข้าหาคำตอบ (Rate of convergence) ที่อยู่ที่ประมาณ 1.84 ซึ่งสูงกว่าวิธีการ Secant ที่อยู่ที่ประมาณ 1.62 ในขณะที่วิธีการ Newton-Raphson นั้นจะอยู่ที่ 2.0 และจุดเด่นอีกข้อหนึ่งของวิธีการ Müller คือการที่สามารถหารากที่เป็นจำนวนเชิงซ้อนได้
   

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

ทีนี้เราลองมาดูตัวอย่างการใช้การ สมการข้างล่างเป็นตัวอย่างหนึ่งของสมการคำนวณค่าสัมประสิทธิ์ความเสียดทาน (friction coefficient) ของการไหลในท่อ

เมื่อ f คือค่าสัมประสิทธิ์ความเสียดทานและ Re คือ Reynolds number โดยจะทำการหาค่า f ที่ Re = 1000

การคำนวณกระทำโดยใช้โปรแกรม Spredsheet OpenOffice 4.1.7 บนระบบปฏิบัติการ 64 bit การคำนวณเริ่มต้นโดยเดาค่า x₀ = 0.1, x₁ = 0.07 และ x₂ = 0.02 ผลการคำนวณแสดงไว้ในตารางในหน้าถัดไป

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

ทศนิยมลงท้ายด้วยเลข 5 จะปัดขึ้นหรือปัดลง (การทำวิทยานิพนธ์ภาคปฏิบัติ ตอนที่ ๙๕) MO Memoir : Friday 30 November 2561

ในการคำนวณด้วยระบบเลขทศนิยม เป็นเรื่องปรกติที่เราพบว่าตัวเลขที่คำนวณได้นั้นเป็นทศนิยมไม่รู้จบ หรือมีหลักทศนิยมมากเกินกว่าที่เราต้องการใช้หรือจัดเก็บได้ ทำให้เราต้องทำการปัดทศนิยมส่วนเกินตัวที่อยู่ด้านท้ายทิ้ง และสิ่งที่ทำโดยทั่วไปก็คือ ถ้าตัวเลขในตำแหน่งทศนิยมที่อยู่ถัดจากตำแหน่งที่เราต้องการเก็บนั้นมีค่าเป็น 4 หรือน้อยกว่า เราก็จะปัดเศษนั้นทิ้งไป กล่าวคือสมมุติว่าเราต้องการเก็บตัวเลขเป็นทศนิยม 2 ตำแหน่ง ตัวเลข 3.124 และ 3.121 ก็จะถูกปัดเป็น 3.12 ทั้งสองจำนวน และในทางกลับกันถ้าตัวเลขในตำแหน่งทศนิยมที่อยู่ถัดจากตำแหน่งที่เราต้องการเก็บนั้นมีค่าเป็น 6 หรือมากกว่า เราก็จะปัดตัวเลขที่อยู่ข้างหน้านั้นขึ้น 1 หน่วย กล่าวคือสมมุติว่าเราต้องการเก็บตัวเลขเป็นทศนิยม 2 ตำแหน่ง ตัวเลข 3.126 และ 3.129 ก็จะถูกปัดเป็น 3.13 ทั้งสองจำนวน
 

ประเด็นที่เป็นคำถามก็คือ ถ้าทศนิยมตำแหน่งนั้นเป็นเลข 5 เราควรปัดเศษ 5 ทิ้งหรือปัดตัวเลขข้างหน้าขึ้น 1 หน่วย กล่าวคือในกรณีของเลข 3.125 เราควรจะปัดเป็น 3.12 (คือปัดทิ้ง) หรือเป็น 3.13 (คือปัดขึ้นหนึ่งหน่วย) และการทำแต่ละแบบนั้นส่งผลอย่างไรหรือไม่อย่างไร ถ้าเรานำตัวเลขที่ผ่านการปัดเศษนั้นไปทำการคำนวณต่อ
 

ก่อนอื่นเราลองมาดูตัวอย่างง่าย ๆ กันก่อน โดยสมมุติว่าเราต้องการคำนวณค่าผลรวม 
  

(1/2) + (2/2) + (3/2) + ... + (18/2) + (19/2) + (20/2)
 

ว่าได้เท่ากับเท่าไรโดยปัดทศนิยมให้เป็นจำนวนเต็ม ซึ่งเราอาจทำการคำนวณด้วยรูปแบบต่าง ๆ ได้ดังนี้

รูปแบบที่ 1 บวกตัวเลขที่เป็นตัวเศษทุกตัวเข้าด้วยกันก่อน จากนั้นจึงหารด้วย 2 
  

(1 + 2 + 3 + ... + 18 + 19 + 20)/2 = 210/2 = 105

รูปแบบที่ 2 หารตัวเลขแต่ละตัวด้วย 2 ก่อนทำการบวกเข้าด้วยกัน เศษที่เหลือจากการหาร (0.5) ให้ทำการปัดลงอย่างเดียว 
  

(1/2) + (2/2) + (3/2) + ... + (18/2) + (19/2) + (20/2) = 0 + 1 + 1 + ... + 9 + 9 + 10 = 100

รูปแบบที่ 3 หารตัวเลขแต่ละตัวด้วย 2 ก่อนทำการบวกเข้าด้วยกัน เศษที่เหลือจากการหาร (0.5) ให้ทำการปัดขึ้นอย่างเดียว 
  

(1/2) + (2/2) + (3/2) + ... + (18/2) + (19/2) + (20/2) = 1 + 1 + 2 + ... + 9 + 10 + 10 = 110

รูปแบบที่ 4 หารตัวเลขแต่ละตัวด้วย 2 ก่อนทำการบวกเข้าด้วยกัน ในกรณีที่มีเศษเหลือ ถ้าตัวเลขหน้า .5 เป็นเลขคี่ให้ทำการปัดขึ้น (เช่น 1.5 ปัดเป็น 2) ถ้าตัวเลขหน้า .5 เป็นเลขคู่ให้ทำการปัดลง (เช่น 4.5 ปัดเป็น 4)
 

(1/2) + (2/2) + (3/2) + ... + (18/2) + (19/2) + (20/2) = 0 + 1 + 2 + ... + 9 + 10 + 10 = 105

รูปแบบที่ 5 หารตัวเลขแต่ละตัวด้วย 2 ก่อนทำการบวกเข้าด้วยกัน ในกรณีที่มีเศษเหลือ ถ้าตัวเลขหน้า .5 เป็นเลขคี่ให้ทำการปัดลง (เช่น 1.5 ปัดเป็น 1) ถ้าตัวเลขหน้า .5 เป็นเลขคู่ให้ทำการปัดขึ้น (เช่น 4.5 ปัดเป็น 5)
 

(1/2) + (2/2) + (3/2) + ... + (18/2) + (19/2) + (20/2) = 1 + 1 + 1 + ... + 9 + 9 + 10 = 105
 

จะเห็นว่ารูปแบบที่ 1 ที่มีการปัดเศษเพียงครั้งเดียว และรูปแบบที่ 4 และ 5 ที่ให้กรณีเศษ .5 มีทั้งการปัดขึ้นและลง (ขึ้นอยู่กับตัวเลขข้างหน้าว่าเป็นเลขคู่หรือเลขคี่) ต่างให้คำตอบที่เหมือนกัน
 

ส่วนกรณีที่ทำการคำนวณได้เศษ .5 แล้วทำการปัดลงเสมอนั้น ทำให้โอกาสปัดเศษลงมีมากกว่าปัดเศษขึ้น เมื่อคำนวณไปเรื่อย ๆ จะทำให้ได้คำตอบที่มีแนวโน้มต่ำกว่าคำตอบที่ถูกต้องมากขึ้นเรื่อย ๆ ในทางกลับกันพอคำนวณได้เศษ .5 แล้วทำการปัดขึ้นเพียงอย่างเดียว ทำให้โอกาสปัดเศษขึ้นมีมากกว่าปัดเศษลง เมื่อคำนวณไปเรื่อย ๆ จะทำให้ได้คำตอบที่มีแนวโน้มสูงกว่าคำตอบที่ถูกต้องมากขึ้นเรื่อย ๆ
 

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

ค่าคลาดเคลื่อน (error) ที่เกิดขึ้นระหว่างการคำนวณสามารถแบ่งออกได้เป็น 2 แบบคือค่าคลาดเคลื่อนการปัดเศษ (round-off error) และค่าคลาดเคลื่อนตัดปลาย (truncation error)
 

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

ค่าคลาดเคลื่อนการปัดเศษเป็นสิ่งที่ติดมากับเครื่องมือ ไม่สามารถหลีกเลี่ยงได้ ในระหว่างการคำนวณค่าความคลาดเคลื่อนนี้อาจจะหักล้างกันไปเอง คืออยู่ในลักษณะปัดขึ้นหรือลงในจำนวนครั้งที่เท่า ๆ กัน ในกรณีเช่นนี้ถ้ามีการคำนวณเกิดขึ้น N ครั้งและเครื่องคำนวณมีค่าความเที่ยง (machine precision) เท่ากับ e_m ค่าความคลาดเคลื่อนทั้งหมดอาจจะมีอันดับ (order) แค่ กล่าวคือถ้าค่า machine precision ของเครื่องคือ 10^-16 และมีการคำนวณ 10⁶ ครั้ง ความคลาดเคลื่อนจะอยู่ที่ระดับ √10⁶ x 10^-16 = 10^-13 หรือสูญเสียทศนิยม 3 ตำแหน่งสุดท้ายไป แต่ถ้าหากการปัดเศษมีแนวโน้มไปในทิศทางใดทิศทางหนึ่ง ค่าความผิดพลาดทั้งหมดจะมีขนาดถึง N.e_m กล่าวคือถ้าค่า machine precision ของเครื่องคือ 10^-16 และมีการคำนวณ 10⁶ ครั้ง ความคลาดเคลื่อนจะอยู่ที่ระดับ 10⁶ x 10^-16 = 10^-10 หรือสูญเสียทศนิยม 6 ตำแหน่งสุดท้ายไป

นอกจากนี้ในการคำนวณบางแบบเช่นการลบเลขสองจำนวนที่มีขนาดใกล้กันมาก ค่าคลาดเคลื่อนการปัดเศษสามารถส่งผลกระทบที่มีนัยสำคัญต่อคำตอบที่ได้จากการคำนวณได้อย่างกระทันหันในการคำนวณเพียงครั้งเดียว ตัวอย่างประเภทนี้ได้แค่การหารากของสมการกำลังสองจากสูตรที่ใช้กันทั่วไปคือ

หรือ 

ปัญหาจะเกิดขึ้นเมื่อ ac << b²

วิธีการคำนวณที่ถูกต้องกว่าคือ

โดยที่รากทั้งสองของสมการคือ

(หมายเหตุ : sgn(b) คือเครื่องหมายของค่า b กล่าวคือ sgn(b) มีค่าเป็น 1 ถ้า b > 0 มีค่าเป็น -1 ถ้า b < 0 และมีค่าเป็น 0 ถ้า b = 0)

เพื่อให้เห็นภาพจะลองทำการหารากของสมการ (จุดที่กราฟตัดแกน x) ของสมการ y = ax² + bx + c จะ เมื่อ a = 1 และ c = 1 โดยที่ b มีค่าต่าง ๆ กัน การคำนวณโดยใช้ spread sheet ของโปรแกรม OpenOffice 4.1.5 โดย

วิธี (ก) คำนวณโดยใช้สูตร ค่าที่ได้แสดงไว้ในตารางที่ ๑

วิธี (ข) คำนวณโดยใช้สูตร

และ และ ค่าที่ได้แสดงไว้ในตารางที่ ๒

ค่า x₁ และ x₂ คือค่าจุดที่กราฟตัดแกน x ที่คำนวณได้ ส่วนค่า y₁ และ y₂ คือค่า y ที่ได้จากการแทนค่า x₁ และ x₂ เข้าไปในสมการ ถ้าหาก x₁ และ x₂ ไม่มีความคลาดเคลื่อนใด ๆ ค่า y₁ และ y₂ ควรจะมีค่าเป็นศูนย์ทั้งสองค่า

ตารางที่ ๑  ค่าคากของสมการที่คำนวณตามวิธีการ (ก) ที่ค่า b ต่าง กัน 

ตารางที่ ๒  ค่าคากของสมการที่คำนวณตามวิธีการ (ข) ที่ค่า b ต่าง กัน

พึงสังเกตว่าการคำนวณตามวิธีการ (ข) (ตารางที่ ๒) จะให้ค่า y ที่ถูกต้องมากกว่าค่า y ที่ได้จากการคำนวณตามวิธีการ (ก) (ตารางที่ ๑) แสดงว่าการคำนวณตามวิธีการ (ข) มีค่าความคลาดเคลื่อนที่ต่ำกว่า