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