โจทย์ผิดหรือถูกคะ?? MO Memoir : Tuesday 15 July 2568

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

วันนี้ก็เลยจะขอลองแก้โจทย์ข้อนี้ดู เพราะเดือนหน้าก็ต้องสอนเรื่องนี้อยู่แล้ว

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

เริ่มจากจัดสมการให้อยู่ในรุป x = f(x) ก่อน

x^(x+2) = 3^4x + 6(9^(x+1))

(x + 2)ln(x) = ln(3^4x + 6(9^(x+1)))

x = ln(3^4x + 6(9^(x+1)))/ln(x) - 2

เนื่องจากมีพจน์ ln(x) ปรากฎ ดังนั้นค่า x ต้องมีค่ามากกว่าศูนย์

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

ตารางข้างล่างคือผลการคำนวณที่ได้ (คำนวนไปว่า 110 รอบ)

xtry	f(xtry) = xcal
2	11.4166655999
11.4166655999	18.6030505575
18.6030505575	25.9647810956
25.9647810956	33.0353033009
33.0353033009	39.5064398927
39.5064398927	45.2217478116
----------	----------
----------	----------
71.9197829085	71.9197829086
71.9197829086	71.9197829086
71.9197829086	71.9197829087
71.9197829087	71.9197829087
71.9197829087	71.9197829087

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

ln(x^(x+2) - 6(9^(x+1))) - ln(3^4x) = 0

ถ้าแทนค่า x ลงไปทางด้านซ้ายของเครื่องหมาย "=" แล้วได้คำตอบเป็นศูนย์ ค่านั้นก็เป็นค่าที่ถูกต้อง

ทีนี้จะลองหาคำตอบด้วยคำสั่ง fzero ของโปรแกรม GNU Octave ดูบ้าง

format longe         กำหนดให้แสดงตัวเลขอย่างละเอียด

>> y = @(x) x^(x+2) - 6*9^(x+1) - 3^(4*x);      กำหนดฟังก์ชัน

>> root = fzero(y,2)          เริ่มหาจากจุด x = 2

error: fzero: zero point is not bracketed            อัลกอริทึมไม่สามารถหาคำตอบได้

error: called from

fzero at line 377 column 7

>> root = fzero(y,5)            ทดลองใหม่โดยเริ่มหาจากจุด x = 5

root = -2.495000000000000e+03         อัลกอริทึมให้คำตอบที่น่าสงสัย (-2495)

>> root = fzero(y,10)          ทดลองใหม่โดยเริ่มหาจากจุด x = 10

root = 7.191978290886503e+01      ได้คำตอบเดียวกับระเบียบวิธีทำซ้ำแบบสืบเนื่องที่เริ่มจาก x = 2

>> root = fzero(y,15)          ทดลองใหม่โดยเริ่มหาจากจุด x = 15

root = -7.350000000000000e+02     อัลกอริทึมให้คำตอบที่น่าสงสัย (-735)

>> root = fzero(y,20)         ทดลองใหม่โดยเริ่มหาจากจุด x = 20

root = -9.800000000000000e+02     อัลกอริทึมให้คำตอบที่น่าสงสัย (-930)

>> root = fzero(y,25)         ทดลองใหม่อีกโดยเริ่มหาจากจุด x = 25

root = 7.191978290886503e+01     ได้คำตอบเดียวกับระเบียบวิธีทำซ้ำแบบสืบเนื่องที่เริ่มจาก x = 2

>> root = fzero(y,30)         ทดลองใหม่อีกโดยเริ่มหาจากจุด x = 30

root = 7.191978290886503e+01     ได้คำตอบเดียวกับระเบียบวิธีทำซ้ำแบบสืบเนื่องที่เริ่มจาก x = 2

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

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

โจทย์ข้อนี้ถ้าหาคำตอบด้วยระเบียบวิธีทำซ้ำแบบสืบเนื่องที่แสดงมาตอนแรก จะพบว่าถ้าเริ่มการคำนวณด้วยการใช้จุดเริ่มต้น x = 5, 15 หรือ 20 ก็จะได้คำตอบเป็น 71.91978291 เสมอ

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

หมายเหตุ การใช้ฟังก์ชั้น fzero ของโปรแกรม GNU Octave เคยเขียนไว้ในบทความเรื่อง
๑. การคำนวณเชิงตัวเลข(๒๕)ตัวอย่างการแก้ปัญหาสมการพีชคณิตไม่เชิงเส้นด้วยFunction fzero ของGNU Octave (MOMemoir : Tuesday 4 February 2563) และ
๒. การคำนวณเชิงตัวเลข(๒๙)เปรียบเทียบการแก้ปัญหาสมการพีชคณิตไม่เชิงเส้นด้วยsolver ของGNU Octave (MOMemoir : Tuesday 3 March 2563)

ซึ่งได้แสดงให้เห็นถึงปัญหาการหาคำตอบขอบสมการพีชคณิตไม่เชิงเส้น

 

การคำนวณเชิงตัวเลข (๓๙) ข้อพึงระวังในการใช้ฟังก์ชันพหุนามในการประมาณค่าในช่วง (๓) MO Memoir : Thursday 24 August 2566

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

รูปที่ ๑ กราฟอุณหภูมิและความหนาแน่นของ CO₂ ที่ความดัน 200 bar.a เขียนด้วยโปรแกรม spread sheet ของ OpenOffice 4.1.14 โดยเลือกชนิดกราฟเป็น x-y scatter และให้ลากเชื่อมต่อจุดด้วยเส้นเรียบ โดยตั้งคุณสมบัติของเส้นเป็น ค่า Default ของโปรแกรม (Cubic spline ความละเอียด 20) ในการลากเส้นเชื่อมจุด จะเห็นว่าในช่วงอุณหภูมิตั้งแต่ 90-500ºC เส้นไม่ค่อยจะราบเรียบเท่าใดนัก ข้อมูลนำมาจากเอกสาร "Properties of Carbon Dioxide" ของ Ihre Messer Group GmbH

สิ่งสำคัญเมื่อได้ค่า y ณ ตำแหน่งที่ต้องการมาแล้ว เราจะมั่นใจได้อย่างไรว่าค่าที่ได้มานั้นมีความน่าเชื่อถือในระดับไหน เพราะจะว่าไปด้วยการเปลี่ยนวิธีการสร้างกราฟหรือเลือกจุดที่ใช้ในการสร้างฟังก์ชันประมาณค่าในช่วง ก็มีสิทธิ์ที่จะได้ผลออกมาแตกต่างกันมากได้ และนี่คือหัวข้อที่นำมาเล่าในวันนี้ โดยขอยกตัวอย่างกรณีความหนาแน่นของแก๊สคาร์บอนไดออกไซด์ (CO₂) ที่ความดัน 200 bar.a โดยจะคำนวณหาค่าที่อุณหภูมิ 120ºC จากจุดข้อมูลที่มีอยู่คือ 50, 60, 70, 80, 90, 100, 200, 300, 400 และ 500ºC

เริ่มจากการนำเอาข้อมูลมาวาดกราฟด้วยโปรแกรม spread sheet ของ OpenOffice 4.1.14 โดยเลือกกราฟชนิด x-y scatter และให้ลากเส้นเรียบเชื่อมต่อจุดด้วยการใช้ค่า default ของโปรแกรม (แบบที่คนส่วนใหญ่ชอบทำกันโดยอาจไม่รู้ว่ามันปรับแต่งได้) ผลที่ได้แสดงไว้ในรูปที่ ๑ จะเห็นว่าเส้นกราฟที่ได้มีปัญหาตรงช่วงที่ระยะหว่างระหว่างจุดมีการเพิ่มขึ้นมากกระทันหัน (90, 100 และ 200ºC) และในช่วงอุณหภูมิสูงนั้นกราฟมีลักษณะเป็นลูกคลื่น

ทีนี้ด้วยข้อมูลชุดเดิม แต่ปรับชนิดของเส้นที่ใช้ลากเชื่อมต่อจุด ผลที่ได้แสดงไว้ในรูปที่ ๒ ซึ่งจะมองเห็นได้ว่าลักษณะเส้นกราฟนั้นมีความสมเหตุสมผลมากกว่า คือลาดลงอย่างเดียว ความขันของเส้นกราฟมีการเปลี่ยนแปลงในทิศทางเดียว คือตั้งชันในช่วงอุณหภูมิต่ำและลดลงในช่วงอุณหภูมิสูง โดยไม่มีลักษณะเป็นลูกคลื่น แม้ว่าทั้งสองแบบจะให้ค่าความหนาแน่นที่อุณหภูมิ 120ºC ใกล้เคียงกัน (425 กับ 420 kg/m³)

รูปที่ ๒ ข้อมูลชุดเดียวกันกับรูปที่ ๑ แต่ใช้ B-Spline ความละเอียดเท่ากับ 20 ลำดับขั้นพหุนามเท่ากับ 3 จะเห็นว่าเส้นมีความโค้งที่ราบเรียบกว่า

แต่ถ้าต้องการค่าแบบละเอียดหรือนำไปเขียนโปรแกรมเพื่อคำนวณค่าความหนาแน่นที่อุณหภูมิต่าง ๆ ก็ควรต้องหาค่าพารามิเตอร์ของฟังก์ชันที่ใช้ประมาณค่า แต่คำถามก็คือควรใช้จุดใดบ้างในการสร้างฟังก์ชันประมาณค่า รูปที่ ๓ แสดงการสร้างฟังก์ชันพหุนามกำลัง 2 สำหรับคำนวณค่าความหนาแน่นที่อุณหภูมิต่าง ๆ กัน โดยเส้นส้มนั้นใช้ข้อมูลที่อุณหภูมิ 90, 100 และ 200ºC ส่วนเส้นสีเขียวนั้นใช้ข้อมูลที่อุณหภูมิ 100, 200 และ 300ºC ฟังก์ชันของทั้งสองเส้นเป็นดังนี้

เส้นสีส้ม y = 0.0271836x^2 - 10.3908x + 1249.26

เส้นสีเขียว y = 0.0078408x^2 - 4.58794x + 862.408

จะเห็นว่าในช่วงอุณหภูมิ 100-200ºC ค่าที่ได้จากการเลือกใช้จุดที่แตกต่างกันนั้นแตกต่างกันอยู่มาก โดยเฉพาะในกรณีของช่วงอุณหภูมิ 183- 200ºC นั้น เส้นสีส้มจะให้ค่าความหนาแน่นที่ต่ำกว่าที่อุณหภูมิ 200ºC ซึ่งไม่ถูกต้อง ส่วนค่าความหนาแน่นที่อุณหภูมิ 120ºC เส้นสีส้มให้ค่าไว้ที่ประมาณ 394 kg/m³ ในขณะที่เส้นสีเขียวให้ค่าประมาณ 425 kg/m³

ด้วยเหตุนี้เวลาสอนนิสิตในเรื่องการประมาณค่าในช่วง จึงบอกกับนิสิตเสมอว่า การคำนวณนั้นมันไม่ได้จบที่แค่ได้ตัวเลข แต่มันต้องทดสอบด้วยว่าตัวเลขที่ได้มานั้นมันน่าเชื่อถือแค่ไหน การทดสอบง่าย ๆ ทำได้ด้วยการพิจารณาภาพรวมการเปลี่ยนแปลงค่าฟังก์ชัน (แบบรูปที่ ๑ และ ๒) หรือลองเลือกใช้จุด (ตำแหน่งและ/หรือจำนวน) ที่แตกต่างกันในการสร้างฟังก์ชันประมาณค่าในช่วง ทำการคำนวณค่า y ณ ตำแหน่ง x ที่ต้องการ แล้วตรวจสอบว่าค่าที่ได้จากฟังก์ชันประมาณค่าที่แตกต่างกันนั้นให้ค่าที่แตกต่างกันมากหรือไม่

รูปที่ ๓ การคำนวณด้วยการใช้ฟังก์ชันพหุนามกำลัง 2 ในการประมาณค่าความหนาแน่นในช่วงอุณหภูมิ 100-200ºC เส้นสีส้มได้จากการใช้ข้อมูลที่อุณหภูมิ 90, 100 และ 200ºC ในการสร้างฟังก์ชัน เส้นเขียวได้จากการใช้ข้อมูลที่อุณหภูมิ 100, 200 และ 300ºC ในการสร้างฟังก์ชัน

ตัวเลขที่เท่ากันแต่ไม่เท่ากัน MO Memoir : Monday 14 August 2566

 วิธีการคำนวณเลขยกกำลังมีด้วยกันหลายวิธี ในกรณีที่เลขยกกำลังนั้นเป็นจำนวนเต็ม เราก็สามารถคำนวณได้ด้วยการคูณกันโดยตรง อย่างเช่น x⁴ ก็จะเท่ากับ x*x*x*x

แต่วิธีการคูณกันโดยตรงนี้ไม่สามารถใช้ได้ในกรณีที่เลขยกกำลังไม่ใช่จำนวนเต็ม ในกรณีเช่นนี้ถ้าหากภาษาคอมพิวเตอร์ที่ใช้นั้นมีคำสั่งยกกำลัง เราก็สามารถใช้คำสั่งนั้นแทน เช่น x⁴ = x**4 หรือ x⁴ = x^4 ซึ่งขึ้นอยู่กับภาษา

แต่ถ้าภาษานั้นไม่มีคำสั่งยกกำลัง เราก็สามารถคำนวณได้ด้วยการหาค่า log ก่อนแล้วค่อยทำการ antilog กลับ ดังนั้นการคำนวณ x⁴ ก็จะออกมาในรูปแบบ x⁴ = exp(4*ln(x)) (ในกรณีของการใช้ natural log)

ในทางทฤษฎีแล้วทั้ง 3 วิธีจะให้คำตอบที่เหมือนกัน แต่เมื่อนำมาคำนวณบนเครื่องคอมพิวเตอร์ที่ต้องมีการแปลงเลขฐาน 10 ให้กลายเป็นเลขฐาน 2 ก่อนเริ่มการคำนวณ และแปลงเลขฐาน 2 ที่คำนวณได้ออกมาเป็นเลขฐาน 10 เพื่อก่ารแสดงผลให้ผู้ใช้เห็น มันจะให้ผลเหมือนกันไหม นี่จึงเป็นที่มาของบทความในวันนี้

การคำนวณในที่นี้ใช้โปรแกรม spreadsheet ของ OpenOffice เวอร์ชัน 4.1.14 บนระบบปฏิบัติการ 64 bit โดยจะคำนวณค่า x⁴ ใน 4 รูปแบบด้วยกันคือ คูณกันโดยตรง (z1), ใช้คำสั่งยกกำลัง (z3), การหาค่า log ก่อนแล้วค่อย anti log กลับด้วยการใช้ natural log (z2) และ log ฐาน 10 (z4) จากนั้นทำการเปรียบเทียบผลที่คำนวณได้เมื่อเทียบกับการคูณกันโดยตรง (z1) ด้วยการหาค่าผลต่าง (zn - z1 ซึ่งถ้าค่าเท่ากันผลลัพธ์ต้องเป็นศูนย์) และการหาค่าสัดส่วน (zn/z1 และ z1/zn ซึ่งถ้าค่าเท่ากันผลลัพธ์ต้องเป็นหนึ่ง)

การทดลองเริ่มจากใช้ค่า x = 0.100000000000000 ก่อน ผลการคำนวณแสดงไว้ในตารางที่ ๑ จะเห็นว่าทั้ง 4 วิธีนั้นให้ผลลัพธ์ที่เหมือนกัน เห็นได้จากการได้ค่าผลต่างเป็น 0 และค่าสัดส่วนเป็น 1

ตารางที่ ๑ ผลการคำนวณ x⁴ ที่ค่า x = 0.100000000000000 ด้วยโปรแกรม OpenOffice 4.1.14

x	1.00000000000000E-001
z1 = x*x*x*x	1.00000000000000E-004
4*ln(x)	-9.21034037197618E+000
z3 = exp(4*ln(x))	1.00000000000000E-004
z3 – z1	0.00000000000000E+000
z3/z1	1.00000000000000E+000
z1/z3	1.00000000000000E+000
z4 = x^4	1.00000000000000E-004
z4 - z1	0.00000000000000E+000
z4/z1	1.00000000000000E+000
z1/z4	1.00000000000000E+000
4*log(x)	-4.00000000000000E+000
z5 = 10^(4*log(x))	1.00000000000000E-004
z5 - z1	0.00000000000000E+000
z5/z1	1.00000000000000E+000
z1/z5	1.00000000000000E+000

ด้วยโปรแกรมตัวเดิม ทำเพียงแค่ป้อนค่า x ค่าใหม่โดยใช้ x = 0.123456000000000 ผลลัพธ์ที่คำนวณได้แสดงไว้ในตารางที่ ๒ จะเห็นว่าในกรณีของการคำนวณด้วยการหาค่า log ด้วย natural log และทำการ antilog กลับนั้น ให้ค่าสัดส่วน z3/z1 = 1 แต่พอคำนวณค่าสัดส่วน z1/z3 กลับได้ค่าเท่ากับ 10

ตารางที่ ๒ ผลการคำนวณ x⁴ ที่ค่า x = 0.123456000000000 ด้วยโปรแกรม OpenOffice 4.1.14

x	1.23456000000000E-001
z1 = x*x*x*x	2.32299784284559E-004
4*ln(x)	-8.36748184679549E+000
z3 = exp(4*ln(x))	2.32299784284559E-004
z3 – z1	0.00000000000000E+000
z3/z1	1.00000000000000E+000
z1/z3	10.00000000000000E+000
z4 = x^4	2.32299784284559E-004
z4 - z1	0.00000000000000E+000
z4/z1	1.00000000000000E+000
z1/z4	1.00000000000000E+000
4*log(x)	-3.63395119348891E+000
z5 = 10^(4*log(x))	2.32299784284559E-004
z5 - z1	0.00000000000000E+000
z5/z1	1.00000000000000E+000
z1/z5	1.00000000000000E+000

ทีนี้พอเปลี่ยนค่า x เป็น 0.123456789123456 ผลลัพธ์ที่คำนวณได้แสดงไว้ในตารางที่ ๓ จะเห็นว่าในกรณีของการคำนวณด้วยการหาค่า log ด้วย natural log และทำการ antilog กลับนั้น ให้ค่าสัดส่วน z1/z3 = 10 แต่พอคำนวณค่าสัดส่วน z3/z1 ได้ค่าเท่ากับ 1 ซึ่งกลับกับกรณีในตารางที่ ๒

ตารางที่ ๓ ผลการคำนวณ x⁴ ที่ค่า x = 0.123456789123456 ด้วยโปรแกรม OpenOffice 4.1.14

x	1.23456789123456E-001
z1 = x*x*x*x	2.32305723727476E-004
4*ln(x)	-8.36745627911360E+000
z3 = exp(4*ln(x))	2.32305723727476E-004
z3 – z1	0.00000000000000E+000
z3/z1	10.00000000000000E+000
z1/z3	1.00000000000000E+000
z4 = x^4	2.32305723727476E-004
z4 - z1	0.00000000000000E+000
z4/z1	1.00000000000000E+000
z1/z4	1.00000000000000E+000
4*log(x)	-3.63394008958575E+000
z5 = 10^(4*log(x))	2.32305723727476E-004
z5 - z1	0.00000000000000E+000
z5/z1	1.00000000000000E+000
z1/z5	1.00000000000000E+000

ทีนี้ลองใหม่ ด้วยการใช้โปรแกรม Excel 2016เปิดไฟล์เดิม (ที่เป็น spreadsheet ของ OpenOffice) แลัวบันทึกเป็น .xlsx จากนั้นลองป้อนค่า x เดียวกันเข้าไป ผลลัพธ์ที่ได้แสดงไว้ในตารางที่ ๔-๖

ตารางที่ ๔ ผลการคำนวณ x⁴ ที่ค่า x = 0.100000000000000 ด้วยโปรแกรม Excel 2016

x	1.00000000000000E-001
z1 = x*x*x*x	1.00000000000000E-004
4*ln(x)	-9.21034037197618E+000
z3 = exp(4*ln(x))	1.00000000000000E-004
z3 – z1	0.00000000000000E+000
z3/z1	1.00000000000000E+000
z1/z3	9.99999999999999E-001
z4 = x^4	1.00000000000000E-004
z4 - z1	0.00000000000000E+000
z4/z1	1.00000000000000E+000
z1/z4	1.00000000000000E+000
4*log(x)	-4.00000000000000E+000
z5 = 10^(4*log(x))	1.00000000000000E-004
z5 - z1	0.00000000000000E+000
z5/z1	1.00000000000000E+000
z1/z5	1.00000000000000E+000

ตารางที่ ๕ ผลการคำนวณ x⁴ ที่ค่า x = 0.123456000000000 ด้วยโปรแกรม Excel 2016

x	1.23456000000000E-001
z1 = x*x*x*x	2.32299784284559E-004
4*ln(x)	-8.36748184679549E+000
z3 = exp(4*ln(x))	2.32299784284559E-004
z3 – z1	0.00000000000000E+000
z3/z1	1.00000000000000E+000
z1/z3	9.99999999999999E-001
z4 = x^4	2.32299784284559E-004
z4 - z1	0.00000000000000E+000
z4/z1	1.00000000000000E+000
z1/z4	1.00000000000000E+000
4*log(x)	-3.63395119348891E+000
z5 = 10^(4*log(x))	2.32299784284559E-004
z5 - z1	0.00000000000000E+000
z5/z1	1.00000000000000E+000
z1/z5	1.00000000000000E+000

 

ตารางที่ ๖ ผลการคำนวณ x⁴ ที่ค่า x = 0.123456789123456 ด้วยโปรแกรม Excel 2016

x	1.23456789123456E-001
z1 = x*x*x*x	2.32305723727476E-004
4*ln(x)	-8.36745627911360E+000
z3 = exp(4*ln(x))	2.32305723727476E-004
z3 – z1	0.00000000000000E+000
z3/z1	9.99999999999999E-001
z1/z3	1.00000000000000E+000
z4 = x^4	2.32305723727476E-004
z4 - z1	0.00000000000000E+000
z4/z1	1.00000000000000E+000
z1/z4	1.00000000000000E+000
4*log(x)	-3.63394008958575E+000
z5 = 10^(4*log(x))	2.32305723727476E-004
z5 - z1	0.00000000000000E+000
z5/z1	9.99999999999999E-001
z1/z5	9.99999999999999E-001

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

การคำนวณเชิงตัวเลข (๓๑) การแก้สมการอนุพันธ์ด้วยฟังก์ชันพหุนาม (๖) MO Memoir : Thursday 2 April 2563

ตรงนี้ขอบันทึกไว้นิดนึง เนื้อหาในบทความชุดนี้นับตั้งแต่ตอนที่ ๑ นั้นอิงมาจากตัวอย่างในหนังสือ "Nonlinear Analysis in Chemical Engineering" ที่เขียนโดย Prof. Bruce A. Finlayson เล่มที่ผมมีนั้นเป็นฉบับพิมพ์เมื่อค.ศ. ๑๙๘๐ เหตุผลที่มีเล่มนี้ก็เพราะต้องใช้ในการทำวิทยานิพนธ์ตอนเรียนปริญญาเอก

  

สมัยเรียนปริญญาตรีได้เรียนเรื่องการแก้ปัญหาสมการอนุพันธ์แยกออกมาเป็นวิชาหนึ่งต่างหาก (๓ หน่วยกิต) และเทคนิคหนึ่งในการแก้ปัญหาที่ได้เรียนก็คือการสมมุติว่าคำตอบนั้นประมาณค่าได้ด้วยอนุกรมอนันต์ (Infinite series) สมัยที่เรียนก็ท่องจำกันแต่ว่าถ้าสมการอนุพันธ์มีหน้าตาแบบนี้ ก็ให้ใช้อนุกรมตัวนี้ในการแก้โจทย์ กว่าจะเข้าใจว่าทำไปต้องเป็นเช่นนั้นก็ผ่านไปหลายปีเหมือนกัน กล่าวคืออนุกรมอนันต์แต่ละอนุกรมนั้นมันมีลักษณะเฉพาะตัวของมัน ถ้าเราเลือกอนุกรมอนันต์ที่มีพฤติกรรมใกล้เคียงกับพฤติกรรมของคำตอบของโจทย์ของเรา การลู่เข้าหาคำตอบก็จะรวดเร็วขึ้น กล่าวคือการได้มาซึ่งตัวเลขค่า y หลังจากแทนค่า x เข้าไปในอนุกรมอนันต์ที่เป็นตำตอบนั้นนั้นไม่จำเป็นต้องใช้จำนวนพจน์มาก ลองดูตัวอย่างนี้ได้ใน Memoir ปีที่ ๗ ฉบับที่ ๘๘๗ วันอาทิตย์ที่ ๙ พฤศจิกายน ๒๕๕๗ เรื่อง "ฟังก์ชันแกมมา (Gamma function) และฟังก์ชันเบสเซล (Bessel function)"

  

คำตอบของสมการที่ใช้ในงานทางวิศวกรรมเคมีนั้นจำนวนไม่น้อยเป็นคำตอบที่มีความสมมาตร กล่าวคือ ณ ตำแหน่ง x ใด ๆ ที่ห่างออกไปจากตำแหน่งจุดสมมาตรไม่ว่าจะเป็นทางซ้ายหรือขวาเป็นระยะที่เท่ากัน จะมีค่า y เท่ากัน ตัวอย่างของโจทย์ประเภทนี้ได้แก่การเปลี่ยนแปลงอุณหภูมิและความเข้มข้นในแนวรัศมีของ tubular reactor (ในกรณีของการทำงานแบบ nonisothermal nonadiabatic) ที่มีความสมมาตรรอบแนวแกนกลางของ reactor การเปลี่ยนแปลงความเข้มข้นของสารในอนุภาคตัวเร่งปฏิกิริยาที่เป็นทรงกลมที่มีความสมมาตรรอบจุดศูนย์กลางของทรงกลม เป็นต้น ดังนั้นจะเป็นการดีกว่าถ้าหากเราใช้ฟังก์ชันที่มีความสมมาตรมาเป็นฟังก์ชันที่ใช้ประมาณค่าคำตอบ

  

ฟังก์ชันพหุนามนั้น ถ้ากำลังสูงสุดเป็นเลขคี่ (odd number) ปลายทั้งสองข้าง (คือเมื่อ x เข้าหา ∞ หรือ -∞ ) จะชี้ไปคนละทางกัน แต่ถ้ากำลังสูงสุดเป็นเลขคู่ (even number) ปลายทั้งสองข้างจะชี้ไปคนละทางกัน (คือข้างหนึ่ง y จะเข้าหา ∞ ในขณะที่อีกข้างหนึ่งจะเข้าหา -∞ ) และถ้าเลขยกกำลังนั้นมีแต่เลขคู่ ฟังก์ชันนั้นก็จะมีความสมมาตรรอบแกนสมมาตร (คือห่างจากตำแหน่ง x ที่เป็นจุดสมมาตรไปทางซ้ายและขวาเป็นระยะ x ที่เท่ากัน จะมีค่า y เท่ากัน) และเราก็สามารถใช้ประโยชน์ตรงจุดนี้ในการนำมาสร้าง orthogonal function ที่มีความสมมาตรเพื่อใช้หาคำตอบของปัญหาที่มีความสมมาตร

  

ในหนังสือของ Finlayson หน้า ๙๔ ได้ให้สมการสำหรับคำนวณหา orthogonal function ในช่วง [0, 1] สำหรับฟังก์ชันพนุนามที่เป็นฟังก์ชันคู่ไว้ดังนี้

ดาวน์โหลดบทความ pdf ได้ที่ลิงก์นี้

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

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

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

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

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