การคำนวณเชิงตัวเลข (๒๖) การแก้ปัญหาสมการอนุพันธ์สามัญ ด้วย 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 ผลการคำนวณจะออกมาเป็นอย่างไร