การคำนวณเชิงตัวเลข (๒๙) เปรียบเทียบการแก้ปัญหาสมการพีชคณิตไม่เชิงเส้นด้วย solver ของ GNU Octave MO Memoir : Tuesday 3 March 2563

ในปีค.ศ. ๑๙๗๖ (พ.ศ. ๒๕๑๙) Ding-Yu Peng และ Donal B. Robinson ได้นำเสนอสมการสภาวะ (equation of state หรือ EOS) ที่ได้มาจากการปรับแต่งสมการสภาวะของ Redlich-Kwong ด้วยรูปสมการดังนี้

โดยที่ P คือความดัน T คืออุณหภูมิ R คือค่าคงที่ของแก๊ส V คือปริมาตรต่อโมล ค่า a และ b เป็นค่าที่ขึ้นอยู่กับอุณหภูมิดังนี้ และมีค่าดังนี้

เมื่อ P_c และ T_c คือค่าที่จุดวิกฤต T_r คือค่า reduced temperature (T/T_c) และ ω คือ acentric factor

การคำนวณค่าปริมาตรจำเพาะของส่วนที่เป็นแก๊สและของเหลวเริ่มด้วยการจัดรูปแบบสมการที่ (1) ใหม่ให้อยู่ในรูป

สำหรับแก๊สเอทิลีน การคำนวณใช้ค่าพารามิเตอร์ดังต่อไปนี้

ก่อนอื่น เราลองมาพิจารณารูปร่างหน้าตาของกราฟของสมการที่ (6) ที่ค่า T และ P ต่าง ๆ ดังแสดงไว้ในรูปที่ ๑
  

ในรูปที่ ๑ เส้นที่ (1) เป็นค่าที่ T = 270 K, P = 37.6274 atm หรือจุดความดันไออิ่มตัว ที่ค่านี้สมการที่ (6) จะมีค่า V ที่ทำให้สมการเป็นจริง (ค่า V ที่ทำให้ค่าฟังก์ชันเป็นศูนย์) ได้ถึง 3 คำตอบ แต่ค่าที่มีความหมายมีเพียง 2 ค่าเท่านั้นคือค่าที่น้อยที่สุดที่เป็นค่าปริมาตรส่วนที่เป็นของเหลว (อยู่ระหว่าง 0.08-0.09) และค่าที่มากที่สุดที่เป็นค่าปริมาณส่วนที่เป็นไอ (ที่ประมาณ 0.32) ส่วนค่าที่อยู่ระหว่างกลางสองค่านั้น (ที่ประมาณ 0.15) ไม่มีความหมายใด ๆ
  

รูปที่ ๑ กราฟของสมการที่ (6) เมื่อ (1) T = 270 K, P = 37.6274 atm (2) T = 274 K, P = 37.6274 atm และ (3) T = 282.4 K, P = 49.938 atm

เส้นที่ (2) เป็นค่าที่ T = 274 K, P = 37.6274 atm หรือที่ค่าอุณหภูมิสูงกว่าอุณหภูมิจุดเดือดที่ความดัน 3.8126 MPa ที่ค่านี้สมการที่ (6) จะมีค่า V ที่ทำให้สมการเป็นจริงเพียงค่าเดียว ซึ่งเป็นค่าของส่วนที่เป็นแก๊ส (ในที่นี้คือ 0.34671) กราฟตรงช่วง V ระหว่าง 0.10-0.11 นั้นไม่ตัดแกน x เพียงแค่วิ่งเข้าหาแล้ววกกลับ
  

ส่วนเส้นที่ (3) เป็นค่าที่ T = 282.4 K, P = 49.938 atm หรือค่าที่จุดวิกฤต (critical point) ที่จุดนี้สมการที่ (6) จะมีเพียงคำตอบเดียว คือส่วนที่เป็น "ของไหล (fluid)" (ที่ประมาณ 0.13) คือไม่สามารถระบุได้ว่าเป็นเฟส "ของเหลว (liquid)" หรือเฟส "แก๊ส (gas)" เพราะไม่มีการแบ่งเฟสกันอย่างชัดเจน

คำตอบของสมการที่ (6) นั้นขึ้นอยู่กับอุณหภูมิและความดัน ที่ความดันต่ำและอุณหภูมิที่สูงพอ สมการที่ (6) จะมีเพียงคำตอบเดียวคือปริมาตรของแก๊ส ที่ค่าอุณหภูมิและความดันที่จุดวิกฤต สมการที่ (6) ก็จะมีเพียงคำตอบเดียวคือส่วนที่เป็นปริมาตรที่จุดวิกฤต และที่ค่าความดันที่สูงและอุณหภูมิที่ต่ำพอ สมการที่ (6) จะมีค่า V ที่ทำให้ "สมการเป็นจริง" ได้ถึง3 ค่า โดยค่า V ที่มีค่ามากที่สุดจะเป็นคำตอบของปริมาตรส่วนที่เป็นแก๊ส ค่า V ที่มีค่าน้อยที่สุดจะเป็นคำตอบของปริมาตรส่วนที่เป็นของเหลว ส่วนค่า V ที่อยู่ระหว่างค่าที่มากที่สุดและค่าที่น้อยที่สุดนั้น แม้ว่าจะเป็นค่า V ที่ทำให้สมการเป็นจริงก็ตาม แต่ในความเป็นจริงจะไม่มีความหมายใด ๆ (คือมันไม่ใช่ปริมาตรของเฟสใดเลย) ค่านี้ต้องโยนทิ้งไป
   

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

กลุ่มที่สองจะทำการยกกำลัง 2 ฟังก์ชันก่อน จากนั้นจึงหาค่า x ที่ทำให้ [f(x)]² = 0 โดยการใช้อัลกอริทึมหาค่าต่ำสุด การยกกำลังสองจะทำให้ค่าฟังก์ชันที่เดิมเป็นลบนั้นมีค่ากลายเป็นบวก ดังนั้นตำแหน่งที่ f(x) ตัดแกน x จะเป็นคำแหน่งที่ [f(x)]² มีค่าต่ำสุด (คือมีค่าเป็นศูนย์)
  

solver ของ GNU Octave ที่นำมาใช้หาคำตอบของสมการที่ (6) ในวันนี้มี 3 ตัวด้วยกัน ตัวแรกคือ "fzero" ที่ใช้สำหรับหาคำตอบระบบที่มีสมการเดียว ตัวที่สองคือ "fsolve" ที่สามารถใช้หาคำตอบของระบบที่มีมากกว่า 1 สมการได้ และตัวที่สามคือ "fminsearch" ที่ใช้สำหรับหาค่าต่ำสุดของฟังก์ชัน ทั้ง 3 ตัวนี้ใช้อัลกอริทึมในการหาที่แตกต่างกัน
  

รูปที่ ๒ ชุดคำสั่งที่ใช้หาคำตอบด้วย solver "fzero" ที่ใช้สำหรับหาคำตอบระบบที่มีสมการเดียว

   

รูปที่ ๓ ชุดคำสั่งที่ใช้หาคำตอบด้วย solver "fsolve" ที่สามารถใช้หาคำตอบของระบบมากกว่า 1 สมการ

รูปที่ ๒ - ๔ เป็นโปรแกรมที่เขียนขึ้นเพื่อใช้หาคำตอบด้วย solver ทั้งสามตัวมื่อ T = 270 K และ P = 37.6274 atm (กราฟเส้น 1 ในรูปที่ ๑) โดยทำการทดลองเริ่มหาคำสั่งจากจุดเริ่มต้นเดียวกันเพื่อดูว่าคำตอบที่ได้จาก solver แต่ละตัวนั้นเป็นอย่างไร คำตอบที่ได้แสดงไว้ในตารางที่ ๑ จะเห็นว่า solver ทั้งสามตัวนั้นสามารถหาค่า V ที่ทำให้สมการเป็นจริงได้ทั้ง 3 ค่า แต่คำตอบที่ได้นั้นขึ้นอยู่กับจุดเริ่มต้นการค้นหา โดยเฉพาะในกรณีของ fzero นั้น พบว่าเมื่อเริ่มต้นคำนวณด้วยค่า V ที่ต่ำเกินไป คำตอบที่ได้นั้นไม่ถูกต้อง
   

รูปที่ ๔ ชุดคำสั่งที่ใช้หาคำตอบด้วย solver "fminsearch" ที่ใช้สำหรับหาค่าต่ำสุดของฟังก์ชัน

ตารางที่ ๑ คำตอบ (V) ที่ได้จาก solver แต่ละตัวเมื่อเริ่มทำการคำนวณที่จุดเริ่มต้นต่าง ๆ

จุดเริ่มต้น	fzero	fsolve	fminsearch
0.040	3.610179364491940E-02*	8.353249907278510E-02	3.185034179687500E-01
0.050	-8.715743984356010E-02*	8.353249905837580E-02	3.184936523437500E-01
0.060	8.353249907278150E-02	8.353249907278470E-02	3.184838867187500E-01
0.083	8.353249907278140E-02	8.353249907277970E-02	8.353405761718750E-02
0.084	8.353249907278140E-02	8.353249907278430E-02	8.352697753906250E-02
0.090	8.353249907278180E-02	8.353249907278190E-02	3.184594726562500E-01
0.110	3.184960937500000E-01	1.507003289001040E-01	3.184960937500000E-01
0.014	1.507003293905560E-01	1.507003293905570E-01	1.506811523437500E-01
0.016	1.507003293905560E-01	1.507003293905570E-01	1.507226562500000E-01
0.200	1.507003293905560E-01	1.507003293900700E-01	3.184692382812500E-01
0.210	3.184753939970270E-01	1.507003293905550E-01	3.184594726562500E-01
0.400	3.184753939970270E-01	3.184753940016410E-01	1.506713867187500E-01
0.500	3.184753939970270E-01	3.184753939978760E-01	3.184814453125000E-01

(* คำตอบที่ผิด)
  

ในกรณีของ fzero นั้นพบว่า เมื่อเริ่มต้นคำนวณจากจุดเริ่มต้น V = 0.05 solver ให้ค่า V ออกมาติดลบ แต่เมื่อตรวจสอบค่าของฟังกชันและรายงานการทำงานของ solver พบว่า ที่จุดดังกล่าวมี f(V) = 3.622620753710365e+17 และ solver รายงานผลการคำนวณออกมาเป็น -5 ซึ่งหมายความว่าค่าที่รายงานนั้นเป็น singular point (ในกรณีนี้น่าจะเป็นจุดที่ฟังก์ชันไม่ต่อเนื่อง) และเมื่อเปลี่ยนจุดเริ่มต้นการคำนวณเป็น V = 0.04 คราวนี้ solver ให้ค่า V ออกมาเป็นบวกที่ถ้าไม่พิจารณาให้ดีอาจจะเข้าใจผิดว่าเป็นคำตอนที่ถูกได้ เพราะเมื่อตรวจสอบค่าของฟังกชันและรายงานการทำงานของ solver พบว่า ที่จุดดังกล่าวมี f(V) = 8.629535666875701e+16 และ solver รายงานผลการคำนวณออกมาเป็น -5 ซึ่งหมายความว่าค่าที่รายงานนั้นเป็น singular point เช่นกัน
  

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

ในความเป็นจริงนั้นปริมาตรของสารต้องมีค่าเป็น "บวก" เสมอ จะมีค่าเป็นศูนย์หรือติดลบไม่ได้ แต่จากตัวอย่างที่ยกมานี้จะเห็นได้ว่า คำตอบที่ได้จาก solver ที่ยกมาเป็นตัวอย่างทดลองใช้นั้นสามารถให้ค่าที่เป็น "บวก" ที่มีทั้งค่าที่ "ถูก" และค่าที่ "ผิด" และบางจังหวะยังอาจให้ค่าที่เป็น "ลบ" ได้อีก

ในส่วนของค่าที่เป็น "บวก" เองนั้น แม้ว่ามันจะเป็นค่าที่ "มีความหมาย" แต่มันอาจไม่ใช่ค่าที่ "ถูกต้องสำหรับการนำไปใช้งานต่อ" เช่นเราต้องการค่าปริมาตรจำเพาะของแก๊สไปใช้งาน แต่กลับได้ค่าปริมาตรจำเพาะของของเหลวแทน
   

ในกรณีของการแก้ปัญหาด้วยเทคนิคที่อิงบนวิธีการของ Newton นั้น สำหรับโจทย์ข้อนี้เป็นที่ทราบกันว่าถ้าเริ่มต้นคำนวณจากค่า V ต่ำ ๆ (เช่นเริ่มจากประมาณค่า b แต่เป็นค่า b ไม่ได้นะ เพราะมันจะเกิดการหารด้วยศูนย์) ก็จะได้ค่าปริมาตรจำเพาะของของเหลว และถ้าเริ่มคำนวณจากค่าเริ่มต้นโดยอิงจาก V = RT/P ก็จะได้ค่าปริมาตรจำเพาะของแก๊ส แต่ในกรณีของ solver "fminsearch" กลับพบว่าไม่เป็นเช่นนั้น
  

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

ค่าที่ทำให้ solver หยุดการคำนวณ

ค่าที่ทำให้ค่าของฟังก์ชันเป็นจริง

ค่าที่เป็นคำตอบที่ถูกต้องของสมการ

ค่าที่เป็นคำตอบที่ถูกต้องสำหรับการนำไปใช้งานต่อ

เคยตรวจสอบกันบ้างไหมครับ :) :) :)