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

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

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

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

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

การคำนวณเชิงตัวเลข (๒๘) การแก้ปัญหาสมการอนุพันธ์สามัญ ด้วย ODE solvers ของ GNU Octave ตอนที่ ๓ MO Memoir : Saturday 15 February 2563

การเกิดปฏิกิริยาใน fixed-bed catalystic reactor นั้น ถ้า heat of reaction มีค่าต่ำ ก็จะประมาณได้ว่า reactor ตัวนั้นทำงานแบบอุณหภูมิคงที่ (isothermal) และถ้า heat of reaction มีค่าไม่สูงมาก ก็จะออกแบบให้ reactor นั้นทำงานแบบ adiabatic (ไม่มีการระบายความร้อนออกหรือให้ความร้อนแก่ reactor โดยตรง) กล่าวคือปล่อยให้ปฏิกิริยาดำเนินไปข้างหน้าเรื่อย ๆ จนได้ค่า conversion ที่ต้องการ แต่ถ้าหากว่าอุณหภูมิของระบบนั้นสูง/ต่ำเกินไปก่อนจะได้ค่า conversion ที่ต้องการ ก็ต้องมีการนำเอา reactor หลายตัวมาต่ออนุกรมกัน โดยมีเครื่องแลกเปลี่ยนความร้อน (ที่ระบาย/ป้อนความร้อน) อยู่ระหว่าง reactor แต่ละตัว

  

ในกรณีของการทำงานแบบ adiabatic นั้น ปฏิกิริยาดูดความร้อน (endothermic reaction) มีแนวโน้มที่จะหยุดตัวเอง เพราะเมื่อปฏิกิริยาดำเนินไปข้างหน้าเรื่อย ๆ ความเข้มข้นของสารตั้งต้นจะลดต่ำลง และอุณหภูมิของระบบก็จะลดต่ำลง สองปรากฏการณ์นี้ส่งผลให้อัตราการเกิดปฏิกิริยาลดต่ำลง 

    

ในกรณีของปฏิกิริยาคายความร้อน (exothermic reaction) นั้นเมื่อปฏิกิริยาดำเนินไปข้างหน้า ความเข้มข้นสารตั้งต้นที่ลดต่ำลงจะดึงให้อัตราการเกิดปฏิกิริยาลดลง แต่อุณหภูมิที่สูงขึ้นจากความร้อนที่คายออกมาจะทำให้ปฏิกิริยาเกิดเร็วขึ้น สองปรากฏการณ์นี้ให้ผลที่ขัดแย้งกัน ถ้าหากการทำปฏิกิริยานั้นเกิดขึ้นในเครื่องปฏิกรณ์แบบ adiabatic (ไม่มีการระบายความร้อนออก) ผลที่เกิดขึ้นจะเกิดในทิศทางเดียวคืออุณหภูมิของระบบจะเพิ่มสูงขึ้นเรื่อย ๆ จนอาจเหนือการควบคุม (ที่เรียกว่าปฏิกิริยาเกิดการ runaway)

   

เพื่อที่จะ "ลด" โอกาสเกิดเหตุการณ์ดังกล่าว ในกรณีของปฏิกิริยาคายความร้อนสูงจึงมักดำเนินการในรูปแบบที่เรียกว่า non-isothermal non-adiabatic คือมีการระบายความร้อนออกจาก reactor โดยตรง ดังนั้นอุณหภูมิใน reactor จึงขึ้นกับอัตราการคายความร้อนจากปฏิกิริยาและอัตราการระบายความร้อนออก เนื่องจากอัตราการระบายความร้อนออกนั้นขึ้นอยู่กับผลต่างระหว่างอุณหภูมิภายใน reactor กับอุณหภูมิแหล่งรับความร้อน ดังนั้นในช่วงที่ปฏิกิริยาเพิ่งเริ่มเกิด (ด้านขาเข้าของ reactor) อุณหภูมิใน reactor (ที่เพิ่มขึ้นจากความร้อนที่ปฏิกิริยาคายออกมา) ยังไม่สูง ผลต่างอุณหภูมิจึงยังไม่สูง อัตราการระบายความร้อนออกจึงต่ำ ส่งผลให้อุณหภูมิใน reactor เพิ่มสูงขึ้นเรื่อย ๆ เมื่อปฏิกิริยาดำเนินไปข้างหน้า แต่เมื่ออุณหภูมิใน reactor สูงถึงระดับหนึ่งจนทำให้ผลต่างอุณหภูมิสูงมากพอ อุณหภูมิใน reactor จะไม่เพิ่มขึ้นและจะเริ่มลดต่ำลง (ผลจากความเข้มข้นสารตั้งต้นที่ลดต่ำลง)

  

ตัวอย่างของปฏิกิริยาประเภทนี้ได้แก่ปฏิกิริยา partial oxidation ไฮโดรคาร์บอนไปเป็นสารประกอบ oxygenate ที่ใช้อากาศเป็นสารออกซิไดซ์ เนื่องจากออกซิเจนในอากาศมีความว่องไวต่ำในการออกซิไดซ์สารตั้งต้น การเริ่มการเกิดปฏิกิริยาจึงต้องการอุณหภูมิที่สูงมากพอ (ไม่เช่นนั้นปฏิกิริยาจะไม่เกิด หรือเกิดช้ามากจนไม่ได้ conversion ที่ต้องการ) อุณหภูมิของแหล่งรับความร้อนจึงต้องสูงตามไปด้วยเพื่อทำให้สารตั้นต้นที่ป้อนเข้ามานั้นมีอุณหภูมิสูงพอที่จะเกิดปฏิกิริยาได้ แต่ในขณะเดียวกันต้องไม่สูงจนทำให้การดึงความร้อนออกจาก reactor ต่ำเกินไป แหล่งรับความร้อนจากปฏิกิริยาที่ใช้กันจึงเป็นพวก molten salt หรือเกลือหลอมเหลวที่เกิดจากสารผสมในสัดส่วนที่เหมาะสมระหว่างเกลือบางชนิด (ทำให้มีจุดหลอมเหลวต่ำเพื่อง่ายต่อการหลอมให้เป็นของเหลว) ตัวอย่างอุณหภูมิของ molten salt นี้อาจอยู่ที่ระดับประมาณ 250-450ºC ขึ้นอยู่กับปฏิกิริยา สภาวะการทำปฏิกิริยา ความว่องไวของตัวเร่งปฏิกิริยา ความเข้มข้นสารตั้งต้น ฯลฯ

ที่ใช้คำว่า "ลด" ในย่อหน้าข้างบนนั่นก็เพราะว่าแม้ว่าจะมีการดึงความร้อนออกจาก reactor โดยตรงมันก็ยังมีโอกาสที่จะเกิดการ runaway อยู่ถ้าหากว่ามีความร้อนสะสมใน reactor มากเกินไป สิ่งที่ทำให้ความร้อนสะสมใน reactor ได้มีทั้ง ความเข้มข้นสารตั้งต้นที่สูง อัตราการไหลที่ต่ำ อุณหภูมิของแหล่งรับความร้อนที่สูงเกินไป อุณหภูมิขาเข้าที่สูงเกินไป และความว่องไวของตัวเร่งปฏิกิริยาที่สูงเกินไป ดังนั้นช่วงการทำงานของ reactor เหล่านี้จึงต้องอยู่ในช่วงที่เหมาะสมที่สามารถทำให้ได้ค่า conversion ที่สูงมากพอโดยไม่เกิดการ runaway

  

การคาดการณ์การเกิด runaway ทำได้ด้วยการแก้สมการดุลมวลสารและพลังงานของระบบ และทดลองปรับเปลี่ยนค่าพารามิเตอร์ต่าง ๆ เพื่อตรวจสอบดูว่าจะเกิดอะไรขึ้น รูปแบบง่าย ๆ รูปแบบหนึ่งของสมการดุลมวลสารและพลังงานแสดงไว้ข้างล่าง (เรียกว่าแบบจำลอง 1 มิติหรือ 1-dimensional model เพราะสมมุติให้การเปลี่ยนแปลงเกิดขึ้นตามความยาวของ reactor เท่านั้น)

  

ในที่นี้พจน์ C คือความเข้มข้น, T คืออุณหภูมิของ reactor, Tc คืออุณหภูมิของแหล่งรับความร้อน, z คือระยะทางตามความยาวของ reactor, exp(a - b/T) คือ rate constant ที่เปลี่ยนตามอุณหภูมิ, H เป็นผลรวมของความร้อนที่เกิด และ U คือค่าสัมประสิทธิ์การถ่ายเทความร้อนระหว่างภายใน reactor กับแหล่งรับความร้อน สมการดุลมวลสารนั้นประกอบด้วยพจน์ของอัตราการหายไปของสารตั้งต้น ในขณะที่สมการดุลความร้อนนั้นประกอบด้วยความร้อนที่เกิดจากปฏิกิริยาและความร้อนที่มีการดึงออก ระบบสมการทั้งสองเป็นระบบสมการอนุพันธ์สามัญที่มีพจน์ที่เป็น non-linear ร่วมอยู่

   

ในที่นี้การหาคำตอบระบบ 2 สมการข้างต้นจะใช้คำสั่ง lsode, ode45 และ ode15s ของ GNU Octave เปรียบเทียบกันโดยจะหาคำตอบในช่วง z = 0.0 - 0.2 คู่มือของโปรแกรมนั้นไม่ได้ให้รายละเอียดว่าคำสั่งดังกล่าวหาคำตอบด้วยวิธีใด คำสั่ง lsode นั้นให้กำหนดช่วงที่ต้องการหาคำตอบและระบุด้วยว่าต้องการซอยช่วงดังกล่าวออกเป็นช่วงย่อยที่ละเอียดมากน้อยเท่าใด (ยิ่งละเอียดมากก็ยิ่งเข้าใกล้คำตอบที่ถูกต้องมากขึ้น) ในที่นี้ได้ทดลองเลือกแบ่งย่อยเป็น 200 ช่วง (ในโปรแกรมต้องใส่ 201 จุด) ค่าพารามิเตอร์ต่าง ๆ ที่ใช้คือ a = 20, b = 14000, H = 10000, T ขาเข้า = 350ºC, C ขาเข้า = 1.0, และ U = 500 ตัวโปรแกรมที่เขียนขึ้นแสดงไว้ในรูปที่ ๑ ส่วนผลการคำนวณที่ได้ที่ใช้ค่า Tc = 414.0, 414.7 และ 414.8ºC แสดงไว้ในรูปที่ ๒ - ๖ โดยรูปที่ ๒ และ ๓ นั้นได้มาจากคำสั่ง lsode รูปที่ ๔ และ ๕ ได้มาจากคำสั่ง ode45 และรูปที่ ๖ ได้มาจากคำสั่ง ode15s

   

ในกรณีของคำสั่ง lsode และ ode45 นั้นพบว่าสองตัวนี้ให้คำตอบที่เหมือนกันคือการทำงานของ reactor จะเป็นถ้าอุณหภูมิของแหล่งรับความร้อนนั้นไม่เกิน 417.7ºC แต่เมื่อเพิ่มอุณหภูมิของแหล่งรับความร้อนเป็น 417.8ºC พบว่าปฏิกิริยาจะเกิดการ runaway ซึ่งเห็นได้จากผลการคำนวณนั้นให้ค่าอุณหภูมิที่เพิ่มสูงขึ้นกระทันหันโดยให้ค่าอุณหภูมิสูงสุดอยู่ที่หลายพันองศาเซลเซียส การแกว่งของคำตอบที่เกิดขึ้นกับกรณีของ lsode นั้นเกิดจากการที่จำนวนจุดคำนวณนั้นต่ำเกินไป (ใช้เพียงแค่ 200 จุด) แต่ที่ไม่พบในกรณีของ ode45 นั้นเป็นเพราะว่า ode45 ใช้การปรับระยะ step size การคำนวณโดยอัตโนมัติเพื่อให้สอดรับกับอัตราการเปลี่ยนแปลงของคำตอบ แต่สิ่งที่เห็นคือระยะห่างของจุดบริเวณนั้นแคบมาก (อยู่ที่ระดับ 10^-7-10^-6) จำนวนจุดการคำนวณตลอดช่วงอยู่ที่เกือบ 35,000 จุด ทำให้การคำนวณใช้เวลานานมากอย่างเห็นได้ชัด

  

แต่พอเปลี่ยนมาใช้ ode15s ปรากฏว่าคำตอบที่ได้นั้นเปลี่ยนไปเป็นคนละเรื่องเลย คือคำสั่งนี้ไม่สามารถพบตำแหน่งที่ปฏิกิริยาเกิดการ runaway ได้ (ตำแหน่งที่อุณหภูมิและความเข้มข้นเปลี่ยนแปลงอย่างรวดเร็ว) คำตอบที่ได้นั้นเสมือนกับว่าปฏิกิริยาแทบจะไม่เกิดแม้ว่าจะใช้อุณหภูมิ Tc สูงกว่า 415.0ºC ก็ตาม

   

รูปที่ ๑ ชุดคำสั่งที่ใช้ในการคำนวณ ข้างบนเป็นของ lsode ส่วนข้างล่างเป็นของ ode45s

ความแตกต่างของตัวอย่างนี้กับตัวอย่างในตอนที่ ๑ และ ๒ อยู่ตรงที่ในตอนที่ ๑ และ ๒ นั้นการเปลี่ยนแปลงอย่างรวดเร็วหรืออย่างกระทันหันเกิดขึ้นจากจุดเริ่มต้นการคำนวณ แต่ในกรณีนี้เกิดขึ้นระหว่างเส้นการคำนวณไปยังจุดปลายทาง ในกรณีของวิธี multi-step ที่มีการใช้จุดข้อมูลย้อนหลังในการคำนวณ ประกอบกับการเลือกระยะ step size อัตโนมัติ จึงเป็นไปได้ว่าการเปลี่ยนแปลงอย่างช้า ๆ ที่เกิดขึ้นอย่างต่อเนื่องมาก่อนหน้า ทำให้ไม่มีการคาดการณ์ว่าจะมีการเปลี่ยนแปลงอย่างกระทันหันเกิดขึ้นต่อหน้า ทำให้การคำนวณเดินข้ามตำแหน่งดังกล่าวไป 

   

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

   

รูปที่ ๒ กราฟการเปลี่ยนแปลองอุณหภูมิภายใน reactor ที่ได้จากคำสั่ง lsode ที่ค่า Tc = 414.0, 414.7 และ 414.8ºC โดยแบ่งช่วงการคำนวณออกเป็น 200 ช่วงย่อย

   

รูปที่ ๓ กราฟการเปลี่ยนแปลงความเข้มข้นของสารตั้งต้นภายใน reactor ที่ได้จากคำสั่ง lsode ที่ค่า Tc = 414.0, 414.7 และ 414.8ºC โดยแบ่งช่วงการคำนวณออกเป็น 200 ช่วงย่อย

   

รูปที่ ๔ กราฟการเปลี่ยนแปลองอุณหภูมิภายใน reactor ที่ได้จากคำสั่ง ode45 ที่ค่า Tc = 414.0, 414.7 และ 414.8ºC

   

รูปที่ ๕ กราฟการเปลี่ยนแปลงความเข้มข้นของสารตั้งต้นภายใน reactor ที่ได้จากคำสั่ง ode45 ที่ค่า Tc = 414.0, 414.7 และ 414.8ºC

   

รูปที่ ๖ กราฟการเปลี่ยนแปลงอุณหภูมิและความเข้มข้นของสารตั้งต้นภายใน reactor ที่ได้จากคำสั่ง ode15s ที่ค่า Tc = 415.0ºC

การคำนวณเชิงตัวเลข (๒๗) การแก้ปัญหาสมการอนุพันธ์สามัญ ด้วย ODE solvers ของ GNU Octave ตอนที่ ๒ MO Memoir : Thursday 13 February 2563

กลับมาต่อกันจากเมื่อวันอังคาร ทบทวนนิดนึงว่าสมการที่นำมาใช้ทดสอบคือสมการอัตราการเกิดปฏิกิริยาของสาร A เมื่อเวลาใด ๆ ที่มีรูปแบบ โดยที่เมื่อ t = 0, A = 1 สมการนี้มีคำตอบคือ เมื่อ A₀ คือค่า A ที่เวลา t = 0 แต่คราวนี้เพิ่มค่า k เป็น 1000 หรือ ๑๐ เท่าของตอนที่แล้ว
   

การคำนวณยังคงใช้ด้วย solver ของ GNU octave ๔ ตัวด้วยกันคือ ode45, ode23, ode15s และ ode15i และเทคนิค implicit Euler ที่เขียนขึ้นมาเองที่ใช้ Δt = 0.001 คงที่ตลอดการคำนวณ และเนื่องด้วยคำตอบที่ได้นั้นมีการเปลี่ยนแปลงอย่างรวดเร็ว จึงได้เขียนกราฟถึงแค่ t = 0.01 ผลการคำนวณที่ได้แสดงไว้ในรูปที่ ๑ - ๕

   

รูปที่ ๑ คำตอบที่ได้จากการคำนวณด้วย solver ode45 เมื่อใช้ k = 1000

จะเห็นได้ว่าในกรณีนี้ ode23 ยังมีปัญหาเรื่องให้ค่าความเข้มข้นออกมาติดลบเกิดขึ้นเร็ว (แค่ประมาณ t = 0.002 ในขณะที่ ode45 และ ode15s แต่ถ้าพิจารณาความใกล้เคียงคำตอบที่ถูกต้องจะเห็นว่า ode15s ให้ค่าที่ดีกว่า ในกรณีของ ode15i นั้นก็มีปัญหาเรื่องค่าความเข้มข้นที่คำนวณได้นั้นออกมาติดลบ แต่ก็ไปเกิดขึ้นแถว ๆ t = 0.017 ส่วนในกรณีของวิธี Implicit Euler นั้นที่เห็นค่าที่คำนวณได้กับคำตอบที่ถูกต้องนั้นห่างกันมากเป็นเพราะการใช้ Δt ที่กว้าง ซึ่งถ้าลดขนาดของ Δt ลงก็จะได้คำตอบที่ใกล้เคียงกับคำตอบที่ถูกต้องมากขึ้น
    

รูปที่ ๒ คำตอบที่ได้จากการคำนวณด้วย solver ode23 เมื่อใช้ k = 1000

  

รูปที่ ๓ คำตอบที่ได้จากการคำนวณด้วย solver ode15s เมื่อใช้ k = 1000

ในทางทฤษฏีนั้นไม่ว่าจะใช้เทคนิคใด ๆ ในการแก้ปัญหา ถ้าใช้ step size ที่มีขนาดเล็กพอก็จะได้คำตอบที่ถูกต้องเหมือนกันหมด เว้นแต่ว่าจะพบกับโจทย์ที่มีปัญหาเรื่องเสถียรภาพสูงมากจนไม่สามารถใช้วิธีตระกูล explicit ได้ ทำให้ต้องใช้เทคนิคตระกูล implicit เท่านั้น
  

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

  

รูปที่ ๔ คำตอบที่ได้จากการคำนวณด้วย solver ode15i เมื่อใช้ k = 1000

   

รูปที่ ๕ คำตอบที่ได้จากการคำนวณด้วยเทคนิค Implicit Euler เมื่อใช้ k = 1000