วิธี Crank-Nicholson แบบปริยาย (Implicit Scheme)

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

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

ในบทความก่อนหน้านี้เกี่ยวกับ Finite Difference Method ได้แสดงให้เห็นว่าวิธีแบบชัดแจ้ง (Explicit) สำหรับการแก้สมการความร้อนด้วยเชิงตัวเลขนั้นต้องใช้ขนาดก้าวเวลาที่เล็กมากอย่างจำกัด ซึ่งกระตุ้นให้พิจารณาวิธีการใหม่ที่สามารถใช้ก้าวเวลาที่ใหญ่ขึ้นได้ แต่ต้องแลกกับงานคำนวณที่มากขึ้นต่อหนึ่งก้าว วิธีการนี้เรียกว่า Crank-Nicholson Scheme

วิธีชัดแจ้งสำหรับสมการความร้อน

วิธีชัดแจ้งใช้ผลต่างไปข้างหน้า (forward difference) สำหรับอนุพันธ์ตามเวลา และผลต่างแบบกึ่งกลาง (centered second difference) สำหรับอนุพันธ์อันดับสองตามตำแหน่ง:

uin+1​−uin/Δt​​=α - ui+1n​−2uin​+ui−1n​​2/(Δx)2

วิธีปริยาย (Implicit Scheme)

Crank-Nicholson Scheme ปรับปรุงวิธีนี้โดยใช้ค่าเฉลี่ยแบบถ่วงน้ำหนักของอนุพันธ์อันดับสองตามตำแหน่งในช่วงเวลา n และ n+1 คำถามที่เกิดขึ้นคือ ค่าของ uuu ที่เวลา n+1 ยังไม่เป็นที่รู้จัก แล้วจะหามาได้อย่างไร? คำถามนี้จะได้รับคำตอบด้านล่าง

อีกคำถามหนึ่งคือควรแบ่งน้ำหนักของอนุพันธ์อันดับสองทั้งสองช่วงเวลาอย่างไร วิธี Crank-Nicholson ใช้การเฉลี่ยแบบ 50-50 แต่ตัวเลือกอื่นก็สามารถใช้ได้ อย่างไรก็ตามเสถียรภาพของวิธีเกี่ยวข้องกับพารามิเตอร์การถ่วงน้ำหนัก θ ในกรณีทั่วไป:

θ=0⇒Explicit

θ=1⇒Fully Implicit 

θ=0.5⇒Crank-Nicholson (กึ่งกลาง)

รูปแบบของ Crank-Nicholson สำหรับสมการความร้อน 1 มิติ

สมการของ Crank-Nicholson สำหรับสมการความร้อนใน 1 มิติเป็นดังนี้:

Δtuin+1​−uin​​=α21​((Δx)2ui+1n+1​−2uin+1​+ui−1n+1​​+(Δx)2ui+1n​−2uin​+ui−1n​​)

หากกำหนดให้

λ=(Δx)2αΔt​

สมการนี้สามารถจัดรูปใหม่เพื่อแยกพจน์ที่ทราบค่าแล้ว (เวลาปัจจุบัน) และพจน์ที่ไม่ทราบค่า (เวลาถัดไป) ได้เป็น:

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

เนื่องจากด้านซ้ายมีสามพจน์ที่ไม่รู้ค่า (เวลาถัดไป) แสดงว่าในแต่ละก้าวเวลา ต้องแก้สมการเชิงเส้นพร้อมกันหลายสมการ สำหรับแต่ละก้าวเวลา nnn จะมี N−2 สมการ (หากมีจุดตารางทั้งหมด N) ซึ่งสอดคล้องกับจุดภายในของตารางพื้นที่ โดยใช้เงื่อนไขขอบแบบ Dirichlet (ค่าคงที่) ทำให้จุดปลายไม่ต้องแก้ในระบบนี้

ระบบสมการเชิงเส้น (Linear System)

ระบบสมการเชิงเส้นสามารถแทนในรูปของเมทริกซ์เป็น:

Aun+1=Bun

โดยที่:

• A และ B เป็นเมทริกซ์แบบสามแนวทแยง (tridiagonal)
• u^n และ u^n+1 คือเวกเตอร์ค่าความร้อนที่จุดภายในของแท่ง ณ เวลา nnn และ n+1n+1n+1
• iii เริ่มตั้งแต่ 1 ถึง N−2

โครงสร้างของเมทริกซ์ A และ B จะเป็น:

A=​1+λ−2λ​000​−2λ​1+λ⋱⋯⋯​0−2λ​⋱−2λ​0​⋯⋯⋱1+λ−2λ​​000−2λ​1+λ​​

B จะมีโครงสร้างเหมือนกัน แต่แทน 1+λ ด้วย 1−λ1 และ −λ/2 ด้วย λ/2​

การแก้ระบบสมการเชิงเส้น

ความพยายามแรกในการแก้ระบบนี้อาจเป็นการอินเวอร์สเมทริกซ์ A เพื่อหา:

u^n+1=A^−1 Bu^n

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

อ้างอิง : Crank-Nicholson Implicit Scheme

จาก https://www.quantstart.com/articles/Crank-Nicholson-Implicit-Scheme/