การคำนวณเชิงตัวเลข (๑) การแก้สมการอนุพันธ์ด้วยฟังก์ชันพหุนาม MO Memoir : Sunday 20 June 2553

เมื่อพฤหัสที่แล้วได้พบกับนักศึกษาปริญญาเอกจากสถาบันแห่งหนึ่ง ซึ่งบังเอิญผมได้รับเชิญไปเป็นอาจารย์ที่ปรึกษาร่วมของเขา งานที่เขาทำนั้นมันคล้ายกับงานวิทยานิพนธ์ที่ผมทำตอนเรียนปริญญาโทและปริญญาเอก คือเป็นการสร้างแบบจำลองทางคณิตศาสตร์ของเครื่องปฏิกรณ์ชนิดเบดนิ่ง (fixed-bed reactor) ซึ่งตอนนั้นผมใช้วิธี orthogonal collocation ร่วมกับ finite element ในการแก้ปัญหาระบบสมการอนุพันธ์ย่อย

สมัยนั้นเวลาทำ simulation แต่ละครั้งก็ต้องเขียนโปรแกรมภาษา FORTRAN กันเอง ต้องแปลงสมการอนุพันธ์ให้อยู่ในรูประบบสมการพีชคณิตที่เหมาะสมเสียก่อน จากนั้นจึงค่อยแก้ปัญหาระบบสมการพีชคณิต ยังดีที่ตอนนั้นยังพอมี NAG Library ซึ่งเป็นซอร์ฟแวร์รวบรวม code ต่าง ๆ สำหรับแก้ปัญหาระบบสมการพีชคณิต ก็เลยไม่ต้องเขียนโปรแกรมส่วนดังกล่าว แต่ก็ต้องรู้ว่า code ตัวไหนเหมาะกับระบบสมการรูปแบบอย่างไรจะได้เลือกใช้ให้ถูก แต่ถึงกระนั้นตัวโปรแกรมก็ยาวร่วม 1000 บรรทัด

ตอนนั้นการทำ dynamic simulation พึ่งอยู่ในช่วงเริ่มต้นการพัฒนา โปรแกรมจำลองกระบวนการทางวิศวกรรมเคมีที่ขายกันอยู่ในขณะนั้นทำได้เพียงแค่ที่ภาวะคงตัว ทั้งนี้เป็นเพราะข้อจำกัดทางด้านคอมพิวเตอร์โดยเฉพาะในเรื่องของหน่วยความจำ ผมเองต้องเขียนโปรแกรมเพื่อทำ dynamic simulation เลยต้องเรียนรู้เทคนิคการแก้ปัญหาที่จะทำให้แก้ปัญหาได้รวดเร็วและใช้หน่วยความจำให้น้อยที่สุด ซึ่งวิธีที่เลือกใช้ตอนนั้นก็คือ orthogonal collocation ร่วมกับ finite element นั่นเอง

เรื่องที่นักศึกษาผู้นั้นมาปรึกษาผมคือต้องการทำความเข้าใจในเรื่องดังกล่าว (การใช้วิธี orthogonal collocation ร่วมกับ finite element ในการแก้ปัญหาระบบสมการอนุพันธ์ย่อย) แต่ปัญหาก็คือดูเหมือนว่าตำราทางด้านนี้ที่เกี่ยวข้องกับทางวิศวกรรมเคมีแทบไม่มีเลย (แม้แต่ตำราภาษาอังกฤษ) และจากประสบการณ์ของผมเองก็พบว่าตำราดังกล่าวเท่าที่มีอยู่นั้น ถ้าผู้อ่านเองไม่มีพื้นฐานทางด้านการคำนวณเชิงตัวเลขมาก่อน ก็จะอ่านเข้าใจยากหรือไม่เข้าใจเลย (ตอนที่ผมไปเรียนนั้นต้องค่อย ๆ แกะสมการทีละบรรทัดกว่าจะเข้าใจ เพราะคนเขียนเองถือว่าคนที่มาอ่านนั้นควรต้องมีความรู้พื้นฐานทางด้านการคำนวณเชิงตัวเลขอยู่ในระดับหนึ่งแล้ว) และบังเอิญตอนที่เขาจบปริญญาตรีมานั้น สถาบันที่เขาจบมาก็ไม่ได้สอนเรื่องนี้ด้วย (เขาบอกผมอย่างนั้น) ก่อนหน้านี้ผมก็เลยแนะนำเอกสารให้ไปอ่านเพื่อช่วยปูพื้นเรื่องการคำนวณเชิงตัวเลขในการแก้ปัญหาระบบสมการพีชคณิตและสมการอนุพันธ์ (ตอนนี้ก็ไม่รู้เหมือนกันว่าเขาอ่านไปถึงไหนแล้ว) มาถึงวันนั้นผมก็เลยช่วยปูพื้นให้ก่อนเบื้องต้นว่าแนวความคิดของวิธีการดังกล่าวมันมาอย่างไร

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

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

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

สมมุติว่าเรามีการนำความร้อนจากผนังด้านหนึ่ง (x =0) ไปยังผนังอีกด้านหนึ่ง (x = 1) ที่ภาวะคงตัว สมการการนำความร้อนผ่านผนังดังกล่าวมีรูปแบบดังนี้

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

T = a₀ + a₁x + a₂x² สมการ (3)

โดยที่ a₀ a₁ และ a₂ คือค่าสัมประสิทธิ์

เพื่อให้สมการที่ (3) นั้นสอดรับกับเงื่อนไขค่าขอบเขต เราจึงนำเงื่อนไขค่าขอบเขตมาแทนลงในสมการที่ (3) จากเงื่อนไขที่ว่า ณ ตำแหน่ง x = 0 ค่า T = 1 เมื่อแทนลงไปในสมการที่ (3) จะได้

a₀ = 1 สมการ (4)

และที่ตำแหน่ง x = 1 ค่า T = 0 เมื่อแทนลงไปในสมการที่ (3) จะได้

0 = 1 + a₁ + a₂

หรือ a₁ = -(1 + a₂) สมการ (5)

แทนค่าสมการที่ (4) และ (5) ลงไปใน (3) จะได้

T = 1 - (1 + a₂)x + a₂x² สมการ (6)

ซึ่งเมื่อจัดรูปแบบใหม่จะได้

T = (1 - x) + (x² - x)a₂ สมการ (7)

โดยการหาอนุพันธ์สมการที่ (6) จะได้

แทนค่าจากสมการที่ (7) (9) และ (10) ลงไปในสมการที่ (2) จะได้

[1 + (1 - x) + (x² - x)a₂](2a₂ ) + [-1 + (2 x -1)a₂]² = 0 สมการ (11)

ในขณะนี้เราติดค่าสัมประสิทธิ์ที่ไม่รู้ค่าอยู่เพียงตัวเดียวคือ a₂ ซึ่งการหาค่า a₂ ทำได้โดยการเลือกตำแหน่งค่า x มาค่าหนึ่ง (การเลือกตำแหน่งนี้ค่อนข้างเป็นอิสระหรือาจกล่าวได้ว่าเลือกตามใจชอบก็ได้ แต่โดยหลักก็คือควรเลือกให้อยู่ในบริเวณที่คำตอบมีการเปลี่ยนแปลงมาก ส่วนการเลือกตำแหน่งจะส่งผลอย่างไรต่อคำตอบนั้นคอยติดตามชมในตอนต่อไป) แล้วนำค่า x ที่เลือกมานั้นแทนลงไปในสมการที่ (11) ก็จะคำนวณหาค่าสัมประสิทธิ์ a₂ ได้

ในครั้งแรกนี้จะลองเลือกที่ตำแหน่ง x = 0.5 ซึ่งเมื่อแทนค่าลงไปในสมการที่ (11) และจัดรูปแบบใหม่จะได้

-0.5 a₂² + 3a₂ + 1 = 0 สมการ (12)

โดยการแก้สมการที่ (12) และนำค่า a₂ ที่ได้ไปคำนวณหาค่า a₁ จากสมการที่ (5) จะได้

a₁ = -0.683375 และ a₂ = -0.316625

หรือ a₁ = -7.316625 และ a₂ = 6.316625

ซึ่งจะทำให้เรามีคำตอบที่เป็นไปได้สองคำตอบคือ

T = 1 - 0.683375x - 0.316625x² สมการ (13)

หรือ T = 1 - 7.316625 x - 6.316625x² สมการ (14)

ที่นี้เราจะรู้ได้อย่างไรว่าสมการไหนเป็นสมการที่ถูกต้อง เราทำได้โดยคำนวณหาค่าฟลักซ์ความร้อน (heat flux) จากสมการ

โดยการใช้สมการที่ (13) จะได้ค่า heat flux ที่ x = 0 มีค่าเป็น 1.36675 และที่ x = 1 มีค่าเป็น 1.316625 ซึ่งแสดงให้เห็นว่าความร้อนไหลไปในทิศทางเดียวกัน (เครื่องหมายของฟลักซ์ความร้อนที่ตำแหน่งทั้งสองเหมือนกัน) แต่ถ้าใช้สมการที่ (14) จะได้ค่า heat flux ที่ x = 0 มีค่าเป็น 14.663325 และที่ x = 1 มีค่าเป็น -5.316625 ซึ่งให้ทิศทางการไหลของความร้อนที่ตำแหน่ง x = 1 ผิดจากที่ควรเป็น (กล่าวคือมีเครื่องหมายตรงข้ามกับที่ตำแหน่ง x = 0) ดังนั้นจึงสรุปได้ว่าสมการที่ถูกต้องคือสมการที่ (13) ส่วนสมการที่ (14) ก็ให้โยนทิ้งไป

การทดสอบว่าสมการที่ (13) นั้นเป็นคำตอบที่แท้จริงของสมการที่ (2) หรือว่าเป็นเพียงคำตอบโดยประมาณนั้นทำได้โดยการแทนค่าสมการที่ (13) กลับเข้าไปในสมการที่ (2) ซึ่งถ้าสมการที่ (13) เป็นคำตอบที่เป็นจริงแล้ว เราจะพบว่าเมื่อแทนค่าสมการที่ (13) กลับเข้าไปในสมการที่ (2) จะต้องได้ค่าฟังก์ชันเป็น 0 ไม่ว่าที่ตำแหน่ง x ใด ๆ ถ้าหากพบว่าค่าที่ได้นั้นไม่ได้เป็น 0 ทุกตำแหน่ง x ค่านั้นจะกล่าวเป็นส่วนตกค้าง (residual) ซึ่งส่วนตกค้างนี้เป็นตัวบอกว่าการประมาณค่านั้นใกล้เคียงความจริงมากเท่าใด ยิ่งส่วนตกค้างมีค่าน้อยเท่าใดก็แสดงว่าการประมาณค่าใกล้เคียงความจริงมากเท่านั้น ค่าอุณหภูมิ (T) ฟลักซ์ความร้อน และส่วนตกค้าง ณ ตำแหน่ง x ต่าง ๆ กันได้รวบรวมเอาไว้ในตารางที่ 1 และรูปที่ 1 ข้างล่าง

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

ฉบับนี้คงพอแค่นี้ก่อน แล้วฉบับหน้าเราค่อยมาดูกันว่าถ้าเราเลือกตำแหน่ง x เป็นค่าอื่นที่ไม่ใช่ 0.5 จะส่งผลต่อคำตอบที่ได้อย่างไร

รูปผมเองสมัยเรียนปริญญาเอก ตัวเล็ก ๆ ที่อยู่ข้างหน้าคือ micro reactor ส่วนตัวใหญ่ข้างหลังคือ pilot plant

o-Xylene to phthalic anhydride pilot plant.

Department of Chemical Engineering and Chemical Technology,

Imperial College of Science, Technology, and Medicine, University of London.

South Kensington campus, London, U.K.

March 1993.