การแก้สมการการแพร่แบบชัดแจ้ง (Solving the Diffusion Equation Explicitly)

โพสต์นี้เป็นส่วนหนึ่งของบทความชุดเกี่ยวกับ Finite Difference Method (วิธีผลต่างจำกัด) ซึ่งประกอบด้วยหัวข้ออื่น ๆ ได้แก่:

• การประมาณอนุพันธ์ด้วยวิธีผลต่างจำกัด
• การแก้สมการการแพร่แบบชัดแจ้ง
• วิธี Crank-Nicolson แบบปริยาย
• ตัวแก้สมการเมทริกซ์สามแนวทแยงด้วยอัลกอริธึมของโธมัส

ในบทความแรกของชุดนี้ ได้แสดงให้เห็นว่าอนุพันธ์ของฟังก์ชันต่อเนื่องสามารถประมาณได้บนโดเมนที่ไม่ต่อเนื่อง ขั้นตอนถัดไปคือการนำอนุพันธ์เหล่านี้ไปใช้กับสมการเชิงอนุพันธ์ย่อยแบบพาราโบลิก (parabolic PDE) ซึ่งสมการความร้อน (Heat Equation) เป็นตัวอย่างมาตรฐาน และเราจะทำให้สมการนี้ไม่ต่อเนื่อง (discretize)

สมการความร้อนและการทำให้ไม่ต่อเนื่อง

สมการความร้อนแบบหนึ่งมิติ มีรูปแบบ:

∂u/∂t ​= D(∂2u/∂x2)​

โดยที่ u(x,t) คืออุณหภูมิ ณ ตำแหน่ง x และเวลา t, D คือสัมประสิทธิ์การแพร่ (diffusivity)

เราจะประมาณอนุพันธ์:

ใช้ Forward Difference กับอนุพันธ์ตามเวลา:

∂u/∂t ​≈ (uin+1​−uin)​​/Δt

ใช้ Centered Difference อันดับสองกับอนุพันธ์ตามตำแหน่ง:

∂x2∂2u​≈(Δx)2ui+1n​−2uin​+ui−1n​​

แทนค่าลงในสมการ จะได้:

Δtuin+1​−uin​​=D⋅(Δx)2ui+1n​−2uin​+ui−1n​​

จัดรูปใหม่เพื่อหา uin+1​:

uin+1​=uin​+(Δx)2DΔt​(ui+1n​−2uin​+ui−1n​)

ให้ α=(Δx)2DΔt​ จะได้:

uin+1​=uin​+α(ui+1n​−2uin​+ui−1n​)

การคำนวณคำตอบ (Computation of Solution)

สมการนี้สามารถนำไปใช้คำนวณ uin+1​ ได้โดยตรงจากค่าของ uin​, ui−1n​,  และ ui+1n​​​ ที่ทราบอยู่แล้ว จึงเรียกว่าสมการ explicit

หากกำหนดเงื่อนไขเริ่มต้น เช่น ให้:

u0=[0,0,1,0,0,1,0,0,1,0,0]

โดยมี “ยอดคลื่นความร้อน” 3 จุด และเลือกพารามิเตอร์:

Δx=1,Δt=0.001,D=1⇒α=1⋅0.00112=0.001\Delta x = 1,\quad \Delta t = 0.001,\quad D = 1 \Rightarrow \alpha = \frac{1 \cdot 0.001}{1^2} = 0.001Δx=1,Δt=0.001,D=1⇒α=121⋅0.001​=0.001

สามารถใช้สูตร:

Δx=1,Δt=0.001,D=1⇒α=121⋅0.001​=0.001

เพื่อคำนวณค่าต่อเนื่องตามเวลาแบบ time marching

ปัญหาเสถียรภาพ (Stability Issue)

ถ้าเลือก Δt=0.01 จะได้ α=0.01 ซึ่งทำให้คำตอบสั่นอย่างรุนแรง:

• รูปแบบของคลื่นที่ได้จะมีค่าขึ้นลงผิดธรรมชาติ
• ไม่สอดคล้องกับคำตอบจริงตามหลักฟิสิกส์ที่แสดงว่าความร้อนควรค่อย ๆ แพร่กระจาย

เกณฑ์เสถียรภาพ (Stability Criterion)

เพื่อให้วิธี explicit นี้เสถียร ต้องใช้เงื่อนไข:

α=(Δx)2DΔt​≤21​

หากไม่เป็นไปตามนี้ คำตอบจะไม่เสถียรและอาจเป็นคลื่นผิดรูป

ตัวอย่างการคำนวณที่เสถียร

เลือก:

Δx=1,Δt=0.001,D=1⇒α=0.001≤0.5

สูตร:

uin+1​=uin​+0.001(ui+1n​−2uin​+ui−1n​)

ผลลัพธ์จะ:

ไม่มีค่าเป็นลบ

สะท้อนลักษณะการแพร่ของความร้อนได้ถูกต้อง

สรุป

แม้ว่าวิธี explicit นี้จะคำนวณง่าย แต่ข้อเสียสำคัญคือข้อจำกัดของ time step:

• ถ้า Δx เล็กมาก (ตาข่ายละเอียด) จะต้องใช้ Δt ที่เล็กมากตาม
• ส่งผลให้การคำนวณใช้เวลานาน

บทความถัดไปจะแนะนำ วิธีปริยาย เช่น Crank-Nicolson เพื่อหลีกเลี่ยงข้อจำกัดนี้

อ้างอิง : Solving the Diffusion Equation Explicitly

จาก https://www.quantstart.com/articles/Solving-the-Diffusion-Equation-Explicitly/