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

คู่มือการใช้งาน GNU Octave 5.1.0 นั้นกล่าวว่ามีฟังก์ชัน (ที่เรียกว่า solver) สำหรับแก้ปัญหาระบบสมการอนุพันธ์สามัญ (Ordinarty Differential Equations หรือที่ย่อกันว่า ODE) ที่ "compatible" กับของ MATLAB อยู่ด้วยกัน ๔ ตัวคือ ode45, ode23, ode15s และ ode15i (ซึ่งต่างเป็นวิธีการที่ใช้กับการแก้ไขปัญหาค่าเริ่มต้น (initial value problem) โดยวันนี้เราจะลองมาดูกันว่าผลการทำงานของแต่ละวิธีมันเป็นอย่างไร (ด้วยการทดสอบกับโจทย์ง่าย ๆ) แต่ก่อนอื่นเรามาทำความรู้จัก solever แต่ะละตัวก่อน โดยเริ่มจากตัวแรก ode45 และ ode23 ที่ผมคัดรายละเอียดจากคู่มือมาสั้น ๆ ดังแสดงข้างล่าง

[t, y] = ode45 (fun, trange, init)

[t, y] = ode45 (fun, trange, init, ode_opt)

[t, y, te, ye, ie] = ode45 (…)

solution = ode45 (…)

ode45 (…)

Solve a set of non-stiff Ordinary Differential Equations (non-stiff ODEs) with the well known explicit Dormand-Prince method of order 4.

[t, y] = ode23 (fun, trange, init)

[t, y] = ode23 (fun, trange, init, ode_opt)

[t, y, te, ye, ie] = ode23 (…)

solution = ode23 (…)

ode23 (…)

Solve a set of non-stiff Ordinary Differential Equations (non-stiff ODEs) with the well known explicit Bogacki-Shampine method of order 3.

ode45 ใช้วิธี Dormand-Prince method อันดับ 4 ในขณะที่ ode23 ใช้วิธี Bogacki-Shampine อันดับ 3 เป็นตัวแก้ปัญหา ทั้งสองวิธีนี้เป็นวิธีที่พัฒนามาจาก Runge-Kutta อันดับ 4 ที่เป็นที่รู้จักกันทั่วไป รายละเอียดของวิธีการทั้งสองเคยเล่าไว้ใน Memoir ปีที่ ๑๐ ฉบับที่ ๑๔๔๘ วันพุธที่ ๔ ตุลาคม ๒๕๖๐ เรื่อง "การคำนวณเชิงตัวเลข (๑๘) เปรียบเทียบระเบียบวิธี Runge-Kutta" ความแตกต่างของสองวิธีการนี้ก็คือ ode23 จะให้คำตอบที่มีระดับความถูกต้องต่ำกว่า ode45 แต่การคำนวณทำได้รวดเร็วกว่า ตรงนี้ก็ขึ้นอยู่กับว่าผู้ใช้ต้องการคำตอบที่ถูกต้องขนาดไหน (เช่นคำตอบต้องถูกต้องถึงนัยสำคัญตำแหน่งที่ ๑๖ เลยหรือไม่) และความเร็วในการได้คำตอบเป็นเงื่อนไขที่สำคัญสำหรับการคำนวณหรือไม่ (เช่นในงานที่ต้องการนำเอาผลการคำนวณไปใช้ในการควบคุมกระบวนการ)

  

[t, y] = ode15s (fun, trange, y0)

[t, y] = ode15s (fun, trange, y0, ode_opt)

[t, y, te, ye, ie] = ode15s (…)

solution = ode15s (…)

ode15s (…)

Solve a set of stiff Ordinary Differential Equations (ODEs) or stiff semi-explicit index 1 Differential Algebraic Equations (DAEs).

ode15s uses a variable step, variable order BDF (Backward Differentiation Formula) method that ranges from order 1 to 5.

solver ตัวที่สามคือ ode15s ที่มีรายละเอียดว่าใช้กับระบบสมการที่ "Stiff" หรือระบบสมการผสมพีชคณิต-อนุพันธ์ (Differential Algebraic Equations ที่ย่อว่า DAE) การแก้ปัญหาใช้วิธี Backward Differential Formula (ที่เป็นหนึ่งในกลุ่มวิธี multi-step) สำหรับผู้ที่สงสัยว่าสมการที่ "Stiff" นั้นเป็นอย่างไรสามารถอ่านได้ใน Memoir ปีที่ ๑๐ ฉบับที่ ๑๔๓๗ วันอังคารที่ ๑๒ กันยายน ๒๕๖๐ เรื่อง "การคำนวณเชิงตัวเลข (๑๗) เปรียบเทียบการแก้ Stiff equation ด้วยระเบียบวิธี Runge-Kutta และ Adam-Bashforth"

[t, y] = ode15i (fun, trange, y0, yp0)

[t, y] = ode15i (fun, trange, y0, yp0, ode_opt)

[t, y, te, ye, ie] = ode15i (…)

solution = ode15i (…)

ode15i (…)

Solve a set of fully-implicit Ordinary Differential Equations (ODEs) or index 1 Differential Algebraic Equations (DAEs).

ode15i uses a variable step, variable order BDF (Backward Differentiation Formula) method that ranges from order 1 to 5.

solver ตัวสุดท้ายคือ ode15i ที่ให้รายละเอียดว่าเป็นวิธี "fully-implicit" แต่มีการใช้วิธี Backward Differential Formula ร่วมด้วย

  

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

  

สมการที่นำมาใช้ทดสอบคือสมการอัตราการเกิดปฏิกิริยาของสาร A เมื่อเวลาใด ๆ ที่มีรูปแบบ dA/dt = -kA โดยที่เมื่อ t = 0, A = 1 และ k = 100 สมการนี้มีคำตอบคือ A = A0 exp (-kt) เมื่อ A₀ คือค่า A ที่เวลา t = 0

  

เพื่อทำการเปรียบเทียบผลการคำนวณด้วย solver ของ GNU octave จึงได้ทำการเขียนคำสั่งเพื่อหาค่า A ที่เวลา t ใด ๆ ด้วยเทคนิค implicit Euler โดยสมการที่ได้คือ A(i+1) = A(i)/(1+kΔt) โดยการคำนวณนี้ใช้ Δt = 0.001 คงที่ตลอดการคำนวณ

  

รูปที่ ๑ - ๓ คือชุดคำสั่งใช้ในการคำนวณ

  

รูปที่ ๑ ชุดคำสั่งสำหรับการคำนวณด้วย solver ode45, ode23 และ ode15s ทั้งสาม solver มีรูปแบบการเรียกใช้ฟังก์ชันที่เหมือนกัน เปลี่ยนเพียงแค่ชื่อของ solver ตัวที่ต้องการใช้ที่บรรทัดที่ 5 การจัดสมการจัดในรูปแบบ dy/dx = f(x,y) (ตรงบรรทัดที่ 4) ผลการคำนวณที่ได้ถูกบันทึกลงไฟล์ "result4.txt"

  

รูปที่ ๒ ชุดคำสั่งสำหรับการคำนวณด้วย solver ode15i ตัวนี้จะจัดสมการแตกต่างไปจากสามตัวแรก คือจัดในรูปแบบ dy/dx + f(x,y) = 0 และจำเป็นต้องมีการให้ค่า x และ dy/dx ที่เวลาเริ่มต้นด้วย ในกรณีของโจทย์ข้อที่ยกมาเป็นตัวอย่างนี้ค่า x ที่เวลาเริ่มต้นคือ A₀ = 1 และค่า dy/dx ที่เวลาเริ่มต้นคือ -kA₀ = -100 (ในที่นี้ใช้ค่า k₁ = 100) แต่ก็อาจใช้คำสั่ง decicช่วยในการหาค่าเริ่มต้นได้ด้วย สำหรับโจทย์ข้อนี้ทดลองดูแล้วไม่เห็นว่ามันจะต่างกัน ดังนั้นคำสั่งนี้สงสัยว่าน่าจะเหมาะสำหรับกรณีของระบบหลายสมการมากกว่า

คำตอบที่ได้จากการคำนวณแสดงไว้ในรูปที่ ๔ - ๘ (การคำนวณกระทำบนระบบปฏิบัติการ 64 บิต)

  

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

  

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

  

รูปที่ ๓ ชุดคำสั่งสำหรับการคำนวณด้วยเทคนิค Implicit Euler ที่ใช้สำหรับเปรียบเทียบคำตอบ

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

รูปที่ ๖ เป็นคำตอบที่ได้จากการใช้ solver ode15i จะเห็นได้ว่าในกรณีของตัวอย่างนี้คำตอบที่ได้อยู่ในระดับเดียวกันกับเมื่อใช้ ode45 เพียงแต่การใช้ ode15i จะใช้ระยะ Δt ที่แคบกว่า และเมื่อเปลี่ยนไปใช้ ode15i กลับพบพฤติกรรมที่ไม่คาดว่าจะพบคือ คำตอบที่ได้มีการแกว่งขึ้นลงในช่วงท้าย (เห็นได้ชัดจากการที่เครื่องหมายของคำตอบมีการสลับบวกลบ) ซึ่งไม่น่าจะเป็นพฤติกรรมที่เกิดขึ้นกับเทคนิคตระกูล implicit ดังผลที่แสดงในรูปที่ ๙ กราฟในรูปที่ ๙ ที่เห็นไม่แนบชิดกับคำตอบที่ถูกต้องเป็นเพราะ Δt ที่ใช้นั้นกว้างเกินไป ถ้าลด Δt ให้มีขนาดเล็กลงก็จะได้ผลที่ใกล้กันเข้าไปอีก

  

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

  

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

การที่ ode15i ที่ซอร์ฟแวร์อ้างว่าใช้เทคนิค fully implicit ในการแก้ปัญหาแล้วคำตอบมีการแกว่งนี้เป็นเรื่องแปลก ทั้ง ๆ ที่เทคนิคตระกูลนี้มีความสามารถในการควบคุมการสะสมของ error ที่ดี เรื่องนี้คงต้องค่อย ๆ ศึกษากันต่อไป

  

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

  

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

ในตอนหน้าเราจะมาดูกันว่า ถ้าเพิ่มค่า ให้สูงขึ้นไปอีก k ผลการคำนวณจะออกมาเป็นอย่างไร

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

การคำนวณเชิงตัวเลข (๒๓) การแก้ปัญหาสมการเชิงอนุพันธ์สามัญปัญหาเงื่อนไขค่าเริ่มต้นด้วยระเบียบวิธี Bogacki-Shampine และ Predictor-Evaluator-Corrector-Evaluator (PECE) MO Memoir : Monday 27 August 2561

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

ในการแก้ปัญหาระบบสมการอนุพันธ์สามัญ ปัญหาเงื่อนไขค่าเริ่มต้นนั้น คงปฏิเสธไม่ได้ว่าระเบียบวิธีหลักที่ใช้กันส่วนใหญ่ในการแก้ปัญหาคือระเบียบวิธีในกลุ่มรุงเก-คุตตาอันดับ ๔ (4th order Runge-Kutta) ที่พัฒนาโดยใช้ระเบียบวิธีนี้เป็นพื้นฐาน ด้วยเหตุผลหลายประการได้แก่ การที่มันเป็นระเบียบวิธีที่สามารถเริ่มต้นการคำนวณด้วยตนเองได้ (เป็น single step method) เป็นระเบียบวิธีโดยแจ้ง (เป็น explicit method) ทำให้ไม่มีปัญหากับสมการที่ไม่เป็นเชิงเส้น (non-linear) ให้ความถูกต้องของคำตอบสูง และมีค่า step size ของการอินทิเกรตที่ค่อนข้างกว้าง แต่ในกรณีที่ต้องการความรวดเร็วสูงในการคำนวณโดยความถูกต้องของคำตอบเป็นรอง ระเบียบวิธีในกลุ่มนี้ก็ไม่ใช่ตัวเลือกที่ดีนัก หรือในกรณีของสมการที่มีปัญหาเรื่อง stiffness หรือไวสูงต่อการสะสมค่าความคลาดเคลื่อนในการคำนวณ ก็อาจจะใช้การไม่ได้เลย
 

ด้วยเหตุนี้แม้แต่ในโปรแกรมสำเร็จรูปเช่น MATLAB ที่นิยมใช้กันอยู่ในขณะนี้ ก็ยังมีระเบียบวิธีสำหรับการแก้ปัญหาระบบสมการอนุพันธ์สามัญ ปัญหาเงื่อนไขค่าเริ่มต้น ให้เลือกเกือบสิบระเบียบวิธี โดย ODE45 ที่ใช้ระเบียบวิธีที่พัฒนาโดยอิงจากรุงเก-คุตตาอันดับ ๔ เป็นฐาน ก็เป็นเพียงแค่ตัวเลือกหนึ่งเท่านั้นเอง
 

วันนี้จะลองกลับมาพิจารณาตัวอย่างที่เคยยกมาใน Memoir ฉบับวันพุธที่ ๑๕ สิงหาคมและวันอาทิตย์ที่ ๑๙ สิงหาคมที่ผ่านมา โดยจะทดลองใช้ระเบียบวิธีอื่นในการแก้ปัญหาดูบ้าง
 

ดาวน์โหลดไฟล์ pdf ฉบับเต็มได้ที่ลิงค์นี้ 

รู้ทันนักวิจัย (๖) บน simulation ทุกอย่างเป็นไปได้หมด (ภาค ๒) MO Memoir : Thursday 6 July 2560

บทความชุดนี้ไม่ได้เขียนให้นักวิจัยอ่านนะครับ แต่เขียนให้คนที่ต้องทำงานเกี่ยวข้องกับนักวิจัยอ่าน

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

๑. ทำนายค่าได้ดี ดูจากอะไร

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

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

ผมเคยตั้งคำถามให้กับนิสิตระดับบัณฑิตศึกษาของภาค (ที่ทำงานอยู่ในภาคอุตสาหกรรมและมาเรียนต่อในช่วงวันหยุด) ลองพิจารณาว่า ถ้าคุณมีแบบจำลองอยู่ ๒ แบบ แบบจำลองที่หนึ่งนั้นสามารถทำนายผลได้อยู่ในช่วง ±5 % กับแบบจำลองที่สองที่ทำนายผลผิดไป -50 % ตลอด คุณจะเอาแบบจำลองไหนสอบวิทยานิพนธ์ และแบบจำลองไหนไปใช้งานจริง คำตอบที่ได้ก็คือเขาจะเอาแบบจำลองแบบแรกไปใช้ในการสอบวิทยานิพนธ์ ด้วยเหตุผลที่ว่ามันดูดีในแง่วิชาการและความชอบของอาจารย์ เพราะมันประมาณค่าได้ใกล้เคียงกับความเป็นจริง แต่ถ้าจะเอาไปใช้งานจริงจะเลือกใช้แบบจำลองที่สอง เพราะถ้ามันคำนวนค่าออกมาได้เท่าใด ก็คูณสองเข้าไปเลยเพื่อให้ได้ค่าที่เป็นจริง ในขณะที่ถ้าเป็นแบบจำลองแบบแรกนั้นมันมีปัญหาตรงที่เมื่อได้ค่าจากการคำนวณแล้ว มันบอกไม่ได้ว่าควรทำอย่างไรก็ค่านั้นต่อ หักออก 5 % บวกเพิ่ม 5 % หรือไม่ต้องทำอะไรเลย เพื่อให้ได้ค่าที่ถูกต้อง
 

การที่จะรู้ว่าแบบจำลองทางคณิตศาสตร์ที่ตั้งขึ้นมานั้นมีความถูกต้องมากน้อยเพียงใด ต้องใช้การเปรียบเทียบกับผลการทดลองจริง ลองดูตัวอย่างในรูปที่ ๑ ในหน้าถัดไปดูนะครับ เป็นกราฟความเข้มข้นสาร B ในปฏิกิริยา A -> B -> C สามเหลี่ยมคือจุดข้อมูลการทดลอง ส่วนเส้นทึบคือค่าที่ได้จากแบบจำลองทางคณิตศาสตร์ ปฏิกิริยารูปแบบนี้ ในช่วงแรกความเข้มข้นสาร B จะเพิ่มขึ้นเรื่อย ๆ เนื่องจากยังมีสาร A ในปริมาณมาก แต่เมื่อไปถึงจุดหนึ่งความเข้มข้นสาร B จะลดลงเนื่องจากการสลายตัวของ B ไปเป็น C นั้นมากกว่าการเกิดสาร B จากสาร A (เพราะสาร A มันหมดไป) บทความดังกล่าวบอกว่าแบบจำลองทางคณิตศาสตร์ของเขานั้นสามารถทำนายผลการทดลอง "ได้ดี" 

รูปที่ ๑ กราฟเปรียบเทียบการทำนายความเข้มข้นของสาร B (ปฏิกิริยา A -> B -> C) ระหว่างข้อมูลการทดลอง (จุดสามเหลี่ยม) กับแบบจำลองทางคณิตศาสตร์ (เส้นทึบ)

ทำนายค่า "ได้ดี" ใช้อะไรเป็นเกณฑ์พิจารณา สำหรับงานทางด้านการเกิดปฏิกิริยาเคมีนั้นการที่จะสร้างแบบจำลองทางคณิตศาสตร์ที่สามารถทำนายค่าได้ตรงกับผลการทดลองแบบแทบไม่ผิดเพี้ยนไม่ว่าจะเป็นสภาวะใดนั้น ผมว่าเป็นเรื่องที่อย่าไปคาดหวังว่าจะได้ เว้นแต่ว่าจะเป็นระบบการเกิดปฏิกิริยาง่าย ๆ ที่กลไกการเกิดปฏิกิริยาไม่ซับซ้อนอะไร ในตัวผมเองนั้นเวลาพิจารณาว่าแบบจำลองทางคณิตศาสตร์นั้นทำนายผลการทดลองได้ดีแค่ไหน ผมจะพิจารณาอยู่ ๓ ประเด็นด้วยกันคือ
 

(ก) ขนาด ของตัวแปรตาม (เช่นความเข้มข้นของสาร การเปลี่ยนแปลงอุณหภูมิของระบบ)
 

(ข) ตำแหน่ง ของตัวแปรอิสระ (เช่นเวลาการทำปฏิกิริยาหรือความยาวของ tubular ractor ที่ให้ความเข้มข้นของสารที่ต้องการมากที่สุดก่อนที่มันจะสลายตัวต่อไป ตำแหน่งใน fixed-bed ที่เกิดจุด hot spot) และ
 

(ค) พฤติกรรมการเปลี่ยนแปลงของตัวแปรตาม (เช่นความเข้มข้นของสารในรูปที่ ๑) เมื่อตัวแปรอิสระ (เช่นเวลา) เปลี่ยนแปลงไป
 

ถ้าแบบจำลองนั้นสามารถให้ความถูกต้องในทั้ง ๓ ประเด็นข้างต้นได้ ก็ถือได้ว่าเป็นแบบจำลองที่ดีมาก แต่สำหรับหลาย ๆ งานแล้ว (โดยเฉพาะแบบจำลองที่เน้นไปที่การควบคุมกระบวนการ) เราอาจต้องให้ความสำคัญกับบางประเด็นเป็นพิเศษ กรณีหนึ่งที่เคยเจอก็คือปฏิกิริยาคายความร้อนที่เกิดขึ้นใน non-isothermal non-adiabatic fixed-bed catalytic reactor ของโรงงานหนึ่ง เครื่องปฏิกรณ์นี้มีความยาวของเบดเกือบ 3 เมตร ปฏิกิริยาคายความร้อนนั้นเป็นปฏิกิริยาที่เร่งตนเองในช่วงแรก เพราะความร้อนที่คายออกมาจะทำให้อัตราการเกิดปฏิกิริยาสูงขึ้น (แม้ว่าความเข้มข้นสารตั้งต้นจะลดลงก็ตาม) พอปฏิกิริยาเริ่มเกิด อุณหภูมิในเบดจะเพิ่มขึ้นอย่างรวดเร็วจนเกิดจุด hot spot ขึ้นในเบด (เรียกว่ากระโดดขึ้นกว่า 100ºC ได้ในระยะทางเพียงแค่ประมาณ 30 cm เท่านั้น ก่อนที่จะลดต่ำลงอันเป็นผลจากสารตั้งต้นหมดไป สิ่งที่เขาต้องการทราบจากแบบจำลองก็คือควรทำการติดตั้งเทอร์โมคับเปิล ณ ตำแหน่งใดใน fixed-bed เพื่อที่จะได้สามารถควบคุมอุณหภูมิไม่ให้เพิ่มมากเกินไป ดังนั้นเขาต้องการแบบจำลองที่สามารถทำนาย "ตำแหน่ง" ที่จะเกิด hot spot ได้ถูกต้องในทุกสภาวะการผลิต ส่วนความถูกต้องของอุณหภูมิสูงสุดที่เกิดที่แบบจำลองทำนายได้นั้นไม่ใช่เรื่องสำคัญ
 

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

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

อีกเรื่องหนึ่งที่ต้องคำนึงในการพิจารณารูปกราฟก็คือ จำนวนจุดข้อมูลที่ใช้วางกราฟและเทคนิคที่ใช้ในการสร้าง interpolation curve เพราะบางทีจุดข้อมูลมันใช้ได้ แต่ด้วยการเลือกเทคนิคการสร้าง interpolation curve ทำให้ได้เส้นกราฟออกมาผิดเพี้ยนไปได้ ดูตัวอย่างเรื่องนี้ได้ใน Memoir ปีที่ ๙ ฉบับที่ ๑๒๕๒ วันอังคารที่ ๑๑ ตุลาคม ๒๕๕๙ เรื่อง "การลาก smooth line เชื่อมจุด (การทำวิทยานิพนธ์ภาคปฏิบัติ ตอนที่ ๘๔)"

๒. ทำไมต้อง ode45

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

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

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

ในวิชาการคำนวณเชิงตัวเลขที่ผมสอนนิสิต ผมจะบอกกับนิสิตเสมอว่า ถ้าคุณมีระบบหลายสมการหลายตัวแปร (ในงานวิศวกรรมเคมี ถ้าจะมีถึงหลักหลายพันสมการหรือมากกว่านั้นก็ไม่ใช่เรื่องแปลก) ถ้าทุกสมการเป็นสมการพีชคณิต (Algebraic Equations) ทั้งหมด หรือเป็นสมการเชิงอนุพันธสามัญ (Ordinary Differential Equations ที่ย่อว่า ODE) ทั้งหมด การแก้ปัญหามันไม่ยากอะไร แต่ถ้าเป็นระบบสมการผสมระหว่างสมการพีชคณิตกับสมการอนุพันธ์ (ที่เรียกว่า Differential Algebratic Equations ที่ย่อว่า DAE) การแก้ปัญหาจะวุ่นวายมากกว่า 
  

เทคนิคการคำนวณเชิงตัวเลข (Numerical methods) ที่ใช้สำหรับแก้ปัญหาระบบสมการ ODE มีอยู่ด้วยกันหลายวิธี แต่ไม่ว่าจะใช้เทคนิคใดก็ตาม ถ้าใช้ step size การคำนวณที่เล็กมากพอ (ขนาดที่เล็กมากพอจะแตกต่างไปตามเทคนิคที่เลือกใช้) ค่าที่คำนวณได้จะลู่เข้าหาคำตอบเดียวกันที่เป็นค่าที่ถูกต้อง ตรงนี้ต่างจากกรณีของสมการอนุพันธ์ย่อย (Partial Differential Equations ที่ย่อว่า PDE) ที่มันมีเรื่อง consistency หรือความคงเส้นคงวาของคำตอบเข้ามาอีก คือใช้ต่างเทคนิคกันอาจได้คำลู่เข้าหาคำตอบที่แตกต่างกัน ทำให้เกิดคำถามได้ว่าคำตอบไหนเป็นคำตอบที่ถูกต้อง

ตารางที่ ๑ คำแนะนำในการเลือก ODE solver ของโปรแกรม MATLAB

ข้อมูลจาก https://www.mathworks.com/help/matlab/math/choose-an-ode-solver.html

Solver	Problem Type	Accuracy	When to Use
ode45	Nonstiff	Medium	Most of the time. ode45 should be the first solver you try.
ode23		Low	ode23 can be more efficient than ode45 at problems with crude tolerances, or in the presence of moderate stiffness.
ode113		Low to High	ode113 can be more efficient than ode45 at problems with stringent error tolerances, or when the ODE function is expensive to evaluate.
ode15s	Stiff	Low to Medium	Try ode15s when ode45 fails or is inefficient and you suspect that the problem is stiff. Also use ode15s when solving differential algebraic equations (DAEs).
ode23s		Low	ode23s can be more efficient than ode15s at problems with crude error tolerances. It can solve some stiff problems for which ode15s is not effective. ode23s computes the Jacobian in each step, so it is beneficial to provide the Jacobian via odeset to maximize efficiency and accuracy. If there is a mass matrix, it must be constant.
ode23t		Low	Use ode23t if the problem is only moderately stiff and you need a solution without numerical damping. ode23 can solve differential algebraic equations (DAEs).
ode23tb		Low	Like ode23s, the ode23tb solver might be more efficient than ode15s at problems with crude error tolerances.
ode15i	Fully implicit	Low	Use ode15i for fully implicit problems f(t,y,y') = 0 and for differential algebraic equations (DAEs) of index 1.

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

แต่เมื่อราคาหน่วยความจำลดต่ำลง และคอมพิวเตอร์มีราคาถูกลงมากเมื่อเทียบกับสมรรถนะที่ได้ เทคนิคที่มักจะเป็นเทคนิคแรกที่เลือกใช้ (และมักจะเป็นเทคนิคสุดท้าย) ในการแก้ปัญหาระบบสมการอนุพันธ์สามัญคือ 4th Runge-Kutta หรือเทคนิคที่มีการพัฒนาเพิ่มเติมโดยใช้ 4th Runge-Kutta เป็นต้นแบบ ซึ่งแม้แต่ในการแก้สมการด้วยการใช้โปรแกรม MATLAB ก็ยังได้แนะนำให้ลองใช้เทคนิคนี้เป็นเทคนิคแรกในการหาคำตอบ ซึ่งมันก็คือ ODE Solver ที่ชื่อ "ode45"
 

ตารางที่ ๑ ข้างต้นเป็นรายชื่อ ODE Solver ที่โปรแกรม MATLAB มีให้เลือกใช้ จะเห็นนะครับว่ามีให้เลือกใช้ตั้ง ๘ ตัว ไม่ได้มีเพียงแค่ ode45 แต่ก็พอจะแบ่งเทคนิคออกได้เป็น ๓ กลุ่มโดยอาศัยเสถียรภาพของการคำนวณเป็นเกณฑ์คือ 
  

(ก) พวกปัญหาแบบทั่วไปที่ไม่ค่อยมีปัญหาใด ๆ เกี่ยวกับเสถียรภาพการคำนวณ (ที่เขาเขียนว่า Nonstiff เทคนิคการแก้ปัญหาจะเป็นวิธีการพวก explicit method และ single step method) 
  

(ข) พวกที่ค่อนข้างจะมีปัญหาเกี่ยวกับเรื่องเสถียรภาพของการคำนวณ (ที่เขาเขียนว่า Stiff เทคนิคการแก้ปัญหาจะเป็นวิธีการพวก explicit method และ multi step method)) และ
 

(ค) พวกที่มีปัญหาเกี่ยวกับเสถียรภาพของการคำนวณมากเช่นพวกระบบสมการ DAE (ที่เขาเขียนว่า Fully implicit เทคนิคการแก้ปัญหาจะใช้วิธีการพวก implicit method และ single step method ซึ่งจะใช้หน่วยความจำมาก แต่สามารถแก้ปัญหาที่เทคนิคในข้อ (ก) และ (ข) ไม่สามารถแก้ได้)
 

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

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

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

อันที่จริงยังมีอีก ๒ เรื่องที่จะเขียน แต่ฉบับนี้ก็เห็นว่าลากยาวมาตั้ง ๕ หน้าแล้ว วันนี้ก็คงต้องขอพอแค่นี้ก่อนครับ