การหาผลเฉลยของสมการเชิงอนุพันธ์สามัญด้วย Python#

บทความนี้กล่าวถึง การเขียนโค้ด Python โดยใช้โมดูล เช่น SymPy และ NumPy / SciPy เพื่อหาผลเฉลยของสมการเชิงอนุพันธ์สามัญ

Keywords: Python, Ordinary Differential Equations (ODE), ODE Solving, Symbolic Methods, Numerical Methods

---

▷ สมการเชิงอนุพันธ์สามัญ#

สมการเชิงอนุพันธ์ (Differential Equations: DE) เป็นสมการที่ประกอบไปด้วยฟังก์ชัน เช่น $y(t)$ ซึ่งเป็นตัวแปรตาม (Dependent Variable) ที่ขึ้นอยู่กับ $t$ ซึ่งเป็นตัวแปรอิสระ หรือ ตัวแปรต้น (Independent Variable) รวมถึงอนุพันธ์อันดับต่าง ๆ ของฟังก์ชัน $y(t)$ เช่น $y'(t)$, $y''(t)$, $...$ และอันดับของสมการเชิงอนุพันธ์คือ อันดับสูงสุดของอนุพันธ์ของฟังก์ชันที่มีอยู่ในสมการ

หากฟังก์ชันขึ้นอยู่กับตัวแปรอิสระเพียงตัวเดียว สมการเชิงอนุพันธ์นี้จะเรียกว่า สมการเชิงอนุพันธ์สามัญ (Ordinary Differential Equation: ODE) แต่ถ้ามีตัวแปรต้นหลายตัวและมีการใช้อนุพันธ์เทียบกับตัวแปรเหล่านั้นหลายตัวจะเรียกว่า สมการเชิงอนุพันธ์ย่อย (Partial Differential Equation: PDE)

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

แต่ในบทความนี้จะกล่าวถึง การใช้ซอฟต์แวร์และคำสั่งต่าง ๆ ของ Python (ใช้เวอร์ชัน Python 3.8) โดยแบ่งเป็นสองแนวทาง คือ

• การประมวลผลเชิงสัญลักษณ์ (Symbolic Methods) โดยใช้คำสั่งของโมดูล SymPy
• การประมวลผลเชิงตัวเลข (Numerical Methods) เช่น การใช้คำสั่งของโมดูล NumPy และ SciPy

ในส่วนของการแสดงผล มีการใช้โมดูล Matplotlib เพื่อแสดงรูปกราฟของฟังก์ชัน และ สำหรับการเขียนและรันโค้ดเพื่อดูผลลัพธ์ ก็สามารถใช้ Python IDE อย่างเช่น Spyder 3 หรือจะใช้ Jupyterlab ก็ได้เช่นกัน

ดังนั้นถ้าจะทดลองรันโค้ด จะต้องมีการติดตั้ง Python 3 Interpreter ในคอมพิวเตอร์ของผู้ใช้ก่อน (Python Download & Install) และติดตั้งโมดูล เช่น matplotlib, numpy, scipy, ipython, jupyterlab โดยใช้คำสั่ง pip

---

▷ ตัวอย่างที่ 1#

ตัวอย่างแรกเป็นโจทย์สมการเชิงอนุพันธ์อันดับที่หนึ่ง $$\frac{dy}{dt} = -2t\, y(t) \\ y'(t) = -2 t\, y(t) = f(t,y) \\ $$

หากคำนวณด้วยวิธีทางคณิตศาสตร์ ก็มีวิธีการดังนี้ โดยใช้เทคนิคที่เรียกว่า วิธีแยกตัวแปร (Separation of Variables) แล้วอินทิเกรตทั้งสองตัวแปรที่เกี่ยวข้อง

$$\begin{align} \frac{1}{y}dy &= -2 t\, dt \\ \int \frac{1}{y}dy &= -2 \int t\, dt \\ ln(y) &= -t^2 + C_0 \\ y(t) &= e^{-t^2 + C_0} \\ &= e^{C_0}\, e^{-t^2} \\ &= C_1\, e^{-t^2} \qquad (C_1 := e^{C_0}) \\ \end{align}$$

โดยที่ $C_0$ และ $C_1$ เป็นค่าคงที่ใด ๆ และไม่เป็นศูนย์ และค่าของ $C_1$ จะหาได้ ถ้ามีการกำหนดเงื่อนไขค่าเริ่มต้น (Initial Condition) เช่น ให้ $y(0) = y_0$ ดังนั้น $$y(0) = C_1\,e^{0} = y_0 \Rightarrow y(t) = y_0\, e^{-t^2} \\ $$

ตัวอย่างการเขียนโค้ด Python เพื่อหาผลเฉลย $y(t)$ ของ ODE โดยวิธีการประมวลผลเชิงสัญลักษณ์ทางคณิตศาสตร์ เริ่มต้นด้วยการสร้างตัวแปรที่เป็นสัญลักษณ์ทางคณิตศาสตร์โดยใช้คำสั่ง symbols() ของ SymPy ได้แก่ t ซึ่งเป็นตัวแปรอิสระ และฟังก์ชันทางคณิตศาสตร์ Function() สำหรับ $y(t)$

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

คำสั่ง dsolve() จะถูกเรียกใช้ เพื่อหาผลเฉลยซึ่งเป็นคำตอบสำหรับฟังก์ชัน $y(t)$ และสามารถแสดงเป็นข้อความเอาต์พุต ได้โดยใช้คำสั่ง print() หรือจะใช้ display() ของ IPython เพื่อแสดงผลลัพธ์ในทางคณิตศาสตร์ให้ดูง่ายขึ้น

ฟังก์ชันเชิงสัญลักษณ์ที่เป็นผลเฉลย จะถูกแปลงให้เป็นฟังก์ชันที่นำไปใช้กับอาร์เรย์ของ numpy เพื่อคำนวณค่าได้ โดยใช้คำสั่ง lambdify()

from sympy import *
import numpy as np
import matplotlib.pyplot as plt
from IPython.display import display

# define a symbol: use t as an independent variable
t, k = symbols('t k')
# define a function y(t)
y = Function('y')(t)
# compute the first derivative of y(t)
dydt = y.diff(t)
# define the ODE
ode  = Eq( dydt, -k*t*y )
# substitute k with 2
ode  = ode.subs( k, 2 )
# set the initial condtion: y(0) = 1.0
ics  = { y.subs(t,0): 1.0 }
# solve the ODE with an initial condition (the IVP)
sol  = dsolve( ode, ics=ics )
print( sol.lhs, '=', sol.rhs )
display( Eq( sol.lhs, sol.rhs ) )

# convert the the solution into a callable function
f = lambdify( t, sol.rhs )
# create an array of 100 time steps between [0,1]
ts = np.linspace(0,1,100)

plt.clf()
plt.plot( ts, f(ts) )
plt.grid()
plt.xlabel('t')
plt.ylabel('y(t)')
plt.show()

รูป: ตัวอย่างการเขียนโค้ด Python และรันโค้ดด้วย Spyder 3 IDE

รูป: กราฟแสดงผลเฉลยสำหรับฟังก์ชัน $y(t)$ ในช่วง $t \in [0,1]$

 

ในเชิงเปรียบเทียบ การเขียนโค้ด Python โดยวิธีการประมวลผลเชิงตัวเลข จะต้องมีการกำหนดช่วงสำหรับตัวแปร $t$ เช่น กำหนดให้ $t \in [t_0, t_f]$ ($t_0=0.0$ และ $t_f=1.0$) เป็นต้น และจะต้องกำหนดค่าตัวเลขสำหรับเงื่อนไขเริ่มต้นด้วย $y(t_0) = y_0$ ซึ่งจะถูกใช้เป็นจุดเริ่มต้นในการคำนวณของค่าของฟังก์ชันในตำแหน่งของ $t$ ถัดไป

การหาผลเฉลยของสมการเชิงอนุพันธ์ด้วย Python เราสามารถใช้คำสั่ง solve_ivp() ของ scipy.integrate ได้ดังนี้

import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

t0 = 0.0  # start point for t
tf = 1.0  # end point for t
h  = 0.05 # step size
y0 = 1.0  # initial value: y(t0) = y0
k  = 2

# the exact solution found by analytic method
f_exact = lambda t,y0: y0*np.exp(-k/2*ts**2)

# y' = f(t,y) = -k*t*y
f = lambda t,y: -k*t*y

ts = np.arange(t0,tf,h) # array of time steps
sol = solve_ivp( f, [t0,tf], [y0], t_eval=ts, method='RK45' )

plt.clf()
plt.plot( ts, sol.y[0], 'o', ts, f_exact(ts,y0), '-' )
plt.grid()
plt.xlabel( 't' )
plt.ylabel( 'y(t)' )
plt.legend( ['numeric','analytic'] )
plt.show()

การใช้คำสั่ง solve_ivp() จะต้องมีอาร์กิวเมนต์ เช่น ฟังก์ชัน f สำหรับการคำนวณค่าของ $f(t,y)$ การกำหนดค่าสำหรับตัวแปร ts (เป็นอาร์เรย์ของ Time Steps) ในช่วง [t0,tf] ซึ่งมีระยะห่างคงที่ h เท่ากับ 0.05
และจะต้องกำหนดค่าเริ่มต้น เช่น [y0]สำหรับ $y(t_0)$ ให้มีค่าเท่ากับ 1.0 ตามลำดับ และมีการเลือกวิธีการในการหาผลเฉลย method ให้เป็น RK45 ซึ่งหมายถึง Runge-Kutta (4-th order)

ผลจากการใช้คำสั่ง solve_ivp() สามารถนำไปใช้ในการแสดงรูปกราฟจากชุดข้อมูล

รูป: แสดงรูปกราฟเปรียบระหว่างผลเฉลยที่ได้จากสองวิธีการ

 

** การคำนวณด้วยวิธีการของออยเลอร์ (Euler Method) **

วิธีการคำนวณเชิงตัวเลข เพื่อหาค่าประมาณของผลเฉลยของสมการเชิงอนุพันธ์อันดับหนึ่ง มีหลายวิธี แต่วิธีการหนึ่งที่เป็นพื้นฐาน เรียกว่า Euler Method แบ่งเป็นสองกรณีคือ

1. Forward (Explicit) Euler Method
2. Backward (Implicit) Euler Method

อนุพันธ์อันดับหนึ่งของฟังก์ชัน $y(t)$ ก็คือ $y'(t)$ หรือ เขียนค่าประมาณได้ดังนี้

$$y'(t) = \frac{d}{dt}y(t) \approx \frac{y(t+h) - y(t)}{h}, \; h \rightarrow 0 \\ y(t+h) \approx y(t) + h\, y'(t) \\ $$

จากสมการเชิงอนุพันธ์อันดับหนึ่ง $y'(t) = f(t,y)$ ก็สามารถเขียนสูตรการคำนวณได้ดังนี้

$$y(t+h) \approx y(t) + h\, f(t,y) \\ $$

ถ้าให้ $h$ เป็นระยะห่างและมีค่าคงที่ (เรียกว่า step-size) การคำนวณจะเริ่มจากจุดแรก คือ $(t_0, y(t_0))$ ไปยังตำแหน่งถัดไป คือ $(t_0+h, y(t_0+h))$ ซึ่งหมายถึง $(t_1, y(t_1))$ ในทิศทางไปข้างหน้า (ค่าของ $t$ เพิ่มขึ้น และ $t_0 < t_1$) และเว้นระยะห่างเท่ากับ $h$

ดังนั้นการคำนวณด้วยวิธี Forward Euler Method จะเป็นดังนี้

$$\begin{align} y_{i+1} &= y_{i} + h\,y'(t_{i}) \\ &= y_{i} + h\,f(t_{i},y_{i}), \;\; i=0,1,2,... \\ \end{align}$$

โดยเริ่มต้นที่จุด $(t_0,y_0), y_0 = y(0)$ และคำนวณค่าถัดไปของฟังก์ชันในตำแหน่งถัดไป $y_1, y_2, ..., y_i$ สำหรับ $t_1, t_2, ..., t_i$ ตามลำดับ ซึ่งมีระยะห่างเท่ากับ $h$

จากโจทย์ในตัวอย่างที่ 1 ก็สามารถเขียนสูตรการคำนวณตามวิธีของออยเลอร์แบบมองไปข้างหน้าได้ดังนี้ $$\begin{align} y'(t) &= -2\, t\, y(t) = f(t,y) \\ y_{0} &= 1 \\ y_{i+1} &= y_{i} - 2 h\, t_{i}\, y_{i}, \;\; i=0,1,2,... \\ \end{align}$$

 

แต่ถ้าเป็นวิธี Backward Euler Method จะเป็นดังนี้

$$\begin{align} y_{i+1} &= y_{i} + h\,y'(t_{i+1}) \\ &= y_{i} + h\,f(t_{i+1},y_{i+1}), \;\; i=0,1,2,... \\ \end{align}$$

ซึ่งจะต้องมีจัดรูปใหม่หรือการแก้สมการ เพื่อหาค่า $y_{i+1}$ ในการคำนวณ และสำหรับโจทย์ตัวอย่างที่ 1 จะเขียนสูตรการคำนวณได้ดังนี้

$$\begin{align} y'(t) &= -2\, t\, y(t) = f(t,y) \\ y_0 &= 1 \\ y_{i+1} &= y_{i} - 2 h\, t_{i+1}, y_{i+1}, \;\; i=0,1,2,... \\ \end{align}$$ และถ้านำไปจัดรูปใหม่ของสมการ จะได้เป็น $$\begin{align} y_{i+1} &= y_{i} - 2 h\, t_{i+1}, y_{i+1} \\ y_{i+1}(1 + 2 h\, t_{i+1}) &= y_{i} \\ y_{i+1} &= \frac{y_{i}}{(1 + 2 h\, t_{i+1})}, \;\; i=0,1,2,... \\ \end{align}$$

ลองมาดูตัวอย่างการเขียนโค้ด Python สำหรับสองวิธีของออยเลอร์ และเปรียบเทียบผลลัพธ์ที่ได้

import numpy as np
import matplotlib.pyplot as plt

# y'(t) = -2*t*y(t) = f(t,y) 
t0 = 0.0  # start point for t
tf = 1.0  # end point for t
h  = 0.05 # step size
y0 = 1.0  # initial value: y(t0) = y0

f = lambda t,y: -2*t*y
f_exact = lambda t: np.exp(-t**2)
ts = np.arange(t0,tf,h)

# forward Euler
y = np.zeros( len(ts) )
y[0] = y0
for i in range(len(ts)-1):
    y[i+1] = y[i] + h*f(ts[i],y[i])
y_f = y

# backward Euler
y = np.zeros( len(ts) )
y[0] = y0
for i in range(len(ts)-1):
    y[i+1] = y[i] / (1+2*h*ts[i+1])
y_b = y

plt.clf()
plt.plot( ts, y_f, 'o', ts, y_b, 'o', ts, f_exact(ts), '-')
plt.grid()
plt.xlabel( 't' )
plt.ylabel( 'y(t)' )
plt.legend(['forward','backward','analytic'])
plt.show()

รูป: แสดงรูปกราฟเปรียบระหว่างผลเฉลยที่ได้จากสองวิธีการของออยเลอร์ (ใช้ค่า $h=0.05$)

 

---

▷ ตัวอย่างที่ 2#

ตัวอย่างที่สองเป็นสมการเชิงอนุพันธ์อันดับที่สองที่มีการกำหนดค่าเริ่มต้นสองเงื่อนไข

$$y'' + 2y' + y = sin(2 t) = f(t,y,y') \\ y(0) = -1, y'(0) = 1 \\ $$

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

$$y'' + 2y' + y = 0 \qquad \text{(homogeneous)}$$

จากสมการที่เป็นเอกพันธ์ สามารถเขียนสมการลักษณะเฉพาะ (Characteristic Equation) ที่เกี่ยวข้อง ซึ่งเป็นพหุนามโดยใช้ตัวแปร $s$ ได้ดังนี้

$$(s^2 + 2s + 1) = (s+1)^2 = 0$$

โดยปรกติแล้ว ถ้าหาราก (roots) ของสมการลักษณะเฉพาะได้เป็น $s_{1,2}$ และไม่เป็นค่าที่ซ้ำกัน ก็สามารถเขียนรูปแบบของฟังก์ชันที่เป็นผลเฉลยได้ดังนี้

$$y_h(t) = C_1\, e^{s_1 t} + C_2\, e^{s_2 t}$$

แต่เนื่องจากในตัวอย่างนี้ ได้รากของสมการที่ซ้ำกัน

$$s_{1,2} = -1$$

ดังนั้นจะต้องเขียนรูปแบบของฟังก์ชัน $y_h(t)$ ที่เป็นผลเฉลยดังนี้ $$y_h(t) = C_1\,e^{s_1 t} + C_2\, t\, e^{s_1 t} = C_1\, e^{-t} + C_2\,t\,e^{-t}$$

ถัดไปเป็นการพิจารณาสมการที่ไม่เป็นเอกพันธ์ และทางขวามือของสมการเป็นฟังก์ชัน $sin(2 t)$ ดังนั้นรูปแบบของฟังก์ชันสำหรับผลเฉลยเฉพาะ เป็นดังนี้

$$y_p(t) = A\, cos(2 t) + B\, sin(2 t) \\ $$

โดยที่ $A$ และ $B$ เป็นค่าคงที่ใด ๆ ดังนั้นจะต้องหาค่า $A$ และ$B$ ที่ทำให้สมการเชิงอนุพันธ์เป็นจริง หากใช้ $y_p(t)$ เป็นผลเฉลยในสมการดังกล่าว

$$\begin{align} y_p(t) &= A\, cos(2 t) + B\, sin(2 t) \\ y'_p(t) &= -2A\, sin(2 t) + 2B\, cos(2 t) \\ y''_p(t) &= -4A\, cos(2 t) - 4B\, sin(2 t) \end{align}$$

จากนั้นให้นำฟังก์ชัน $y_p(t)$ และอนุพันธ์ของฟังก์ชัน $y'_p(t)$ และ $y''_p(t)$ ไปแทนลงในสมการเชิงอนุพันธ์ แล้วจัดกลุ่มแยกตามพจน์สำหรับ $sin(2t)$ กับ $cos(2t)$ ดังนี้

$$(-4A -3B)\,sin(2 t) + (-3A + 4B)\,cos(2 t) = sin(2 t) \\ -4A - 3B = 1 \\ -3A + 4B = 0$$

เมื่อแก้ระบบสมการสองตัวแปรจะได้ $A=-4/25$ และ $B=-3/25$

 

ผลเฉลยทั่วไปที่ได้ $y(t)$ เป็นผลรวมของ $y_h(t)$ และ $y_p(t)$ ในรูปแบบของฟังก์ชันต่อไปนี้

$$\begin{align} y(t) &= (C_1 + C_2\, t)\,e^{-t} -\frac{4}{25}cos(2t) - \frac{3}{25}sin(2t) \\ y'(t) &= (-C_1 + C_2 - C_2 t)\,e^{-t} + \frac{8}{25}sin(2t) - \frac{6}{25}cos(2t) \end{align}$$

จากเงื่อนไขค่าเริ่มต้น $y(0) = -1$ และ $y'(0) = 1$ จะได้ระบบสมการที่สามารถนำไปหาค่า ของ $C_1$ และ $C_2$ ได้ $$\begin{align} y(0) &= C_1 - \frac{4}{25} = -1 \\ y'(0) &= -C_1 + C_2 - \frac{6}{25} = 1 \\ \Rightarrow\; C_1 &= -\frac{21}{25}, \, C_2 = \frac{2}{5} \end{align}$$

โดยสรุป จะได้ฟังก์ชัน $y(t)$ ที่เป็นผลเฉลยของปัญหาค่าเริ่มต้นดังนี้ $$y'' + 2y' + y = sin(2 t), \; y(0) = -1, y'(0) = 1 \\ \begin{align} y(t) &= (-\frac{21}{25} + \frac{2}{5} \, t)\,e^{-t} -\frac{3}{25}sin(2t) -\frac{4}{25}cos(2t) \\ &= \frac{1}{25}(10\, t - 21)\,e^{-t} - \frac{1}{5}cos(2t - atan(\frac{3}{4})) \end{align}$$

จะเห็นได้ว่า การหาผลเฉลยของสมการเชิงอนุพันธ์ในตัวอย่างนี้มีหลายขั้นตอน

 

ถัดไปลองมาดูตัวอย่างเขียนโค้ด Python ด้วยวิธีการคำนวณเชิงสัญลักษณ์ เพื่อใช้ในการหาคำตอบ

from sympy import *
import numpy as np
import matplotlib.pyplot as plt
from IPython.display import display, Latex

# define a symbol: use t as an independent variable
t = symbols('t')
# define a function y(t)
y = Function('y')(t)
# compute the first derivate of y(t)
D1y = y.diff(t)
D2y = y.diff(t,2)

# 2nd-order ODE: y'' + 2y' + y = sin(2t)
# Initial conditions: y(0)=-1, y'(0)=1
# define the ODE
ode  = Eq( D2y + 2*D1y + y, sin(2*t) )
# set the initial condition: y(0) = 1.0
ics  = { y.subs(t,0): -1, D1y.subs(t,0): 1 }
# solve the ODE with an initial condition (the IVP)
sol  = dsolve( ode, ics=ics )
print( sol.lhs, '=', sol.rhs )
display( Eq( sol.lhs, sol.rhs ) )
# => y(t) = (2*t/5 - 21/25)*exp(-t) - 3*sin(2*t)/25 - 4*cos(2*t)/25

# convert the the solution into a callable function
f = lambdify( t, sol.rhs )
# create an array of 100 time steps between [t0,tf]
t0 = 0.0
tf = 10.0
ts = np.linspace(t0,tf,100)

plt.clf()
plt.plot( ts, f(ts) )
plt.grid()
plt.xlabel('t')
plt.ylabel('y(t)')
plt.show()

รูป: ตัวอย่างการรันโค้ดใน Spyder IDE สำหรับตัวอย่างที่ 2

รูป: แสดงรูปกราฟสำหรับ $y(t)$ ที่ได้จากการคำนวณเชิงสัญลักษณ์ด้วย SymPy

 

การหาผลเฉลยสำหรับสมการเชิงอนุพันธ์ที่มีอันดับมากกว่าหนึ่ง มีอีกวิธีหนึ่งคือ การแปลงสมการดังกล่าวให้เป็นระบบสมการอนุพันธ์ที่มีอันดับเป็นหนึ่ง ซึ่งมีแนวทางดังนี้

ขั้นตอนแรกคือ กำหนดฟังก์ชันเป็นตัวช่วย เช่น $z_1(t)$ และ $z_2(t)$ สำหรับสมการอนุพันธ์อันดับสอง

$$z_1(t) := y(t), \; z_2(t) := y'(t) \\ z'_1(t) = z_2(t), \; z'_2(t) = y''(t) \\ z_1(0) = y(0), \; z_2(0) = y'(0) \\ $$

จากนั้นเขียนสมการอนุพันธ์ให้อยู่ในรูปเมทริกซ์

$$\begin{align} y'' + 2y' + y &= sin(2t) \\ y'' &= -2y' - y + sin(2t) = f(t,y,y') \\ &= -2 z_2(t) - z_1(t) + sin(2t) \end{align}\\ $$

$$\frac{d}{dt} \begin{pmatrix} z_1(t) \\ z_2(t) \\ \end{pmatrix} = \begin{pmatrix} 0 & 1 \\ -1 & -2 \\ \end{pmatrix} \begin{pmatrix} z_1(t) \\ z_2(t) \\ \end{pmatrix} + \begin{pmatrix} 0 \\ sin(2t) \\ \end{pmatrix}$$

 

ตัวอย่างการเขียนโค้ด Python โดยใช้การคำนวณเชิงสัญลักษณ์ในรูปแบบเมทริกซ์และเวกเตอร์ มีดังนี้

from sympy import *
from sympy.solvers.ode.systems import dsolve_system
import numpy as np
import matplotlib.pyplot as plt
from IPython.display import display

# define a symbol: use t as an independent variable
t = symbols('t')
# define a function y(t)
y, z1, z2 = symbols('y z1 z2', cls=Function)
# compute the first and second derivatives of y(t)
D1y = Derivative(y(t),t,1)
D2y = Derivative(y(t),t,2)

# 2nd-order ODE: y'' + 2y' + y = sin(2t)
# Initial conditions: y(0)=-1, y'(0)=1
# define the ODE (system of 1st-order differential equations)
# [z1'(t), z2'(t)] = [ z2(t), f(t,z1,z2) ]
ode  = [ Eq( z1(t).diff(t), z2(t) ), 
         Eq( z2(t).diff(t), -2*z2(t) -z1(t) + sin(2*t) ) ]
# set the initial condition: y(0) = 1.0
ics  = { z1(0): -1, z2(0): 1 }
# solve the ODE with an initial condition (the IVP)
sol = dsolve( ode, ics=ics )
sol_eq = expand( sol[0].rhs ) 
display( sol_eq )

# convert the the solution into a callable function
f = lambdify( t, sol_eq )
# create an array of 100 time steps between [t0,tf]
t0 = 0.0
tf = 10.0
ts = np.linspace(t0,tf,100)

plt.clf()
plt.plot( ts, f(ts), 'o' )
plt.grid()
plt.xlabel('t')
plt.ylabel('y(t)')
plt.show()

รูป: ผลลัพธ์จากการรันโค้ดเพื่อสาธิตการใช้คำสั่ง dsolve()

 

ลองมาดูตัวอย่างการคำนวณด้วยวิธีเชิงตัวเลข โดยใช้คำสั่ง solve_ivp() ของ SciPy

import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

# 2nd-order ODE: y'' + 2y' + y = sin(2t)
# Initial conditions: y(0)=-1, y'(0)=1
t0 = 0.0   # start point for t
tf = 10.0  # end point for t
h  = 0.05  # step size

f = lambda t,z: [z[1], -2*z[1] - z[0] + np.sin(2*t) ]
ts  = np.arange(t0,tf,h) # array of time steps
sol = solve_ivp( f, [t0,tf], [-1,1], t_eval=ts, method='RK45' )

plt.clf()
plt.plot( ts, sol.y[0], 'o' )
plt.grid()
plt.xlabel( 't' )
plt.ylabel( 'y(t)' )
plt.show()

รูป: แสดงรูปกราฟของฟังก์ชันที่เป็นผลเฉลยเมื่อใช้คำสั่ง solve_ivp()

 

---

▷ ตัวอย่างที่ 3#

ตัวอย่างถัดไปเป็นสมการเชิงอนุพันธ์อันดับ 4 พร้อมเงื่อนไขค่าเริ่มต้น

$$y^{(4)}(t) + 2 y^{(3)}(t) + 2 y''(t) + 2 y'(t) + y(t) = 10 e^{-2t} \\ y(0)=0, y'(0)=1, y''(0)=-1, y'''(0)=2$$

หากพิจารณาสมการที่เป็นเอกพันธ์ จะเขียนสมการลักษณะเฉพาะได้ดังนี้

$$s^4 + 2 s^3 + 2 s^2 + 2 s + 1 = 0 \\ \Rightarrow\; (s+1)^2(s^2 + 1 ) = 0 \\ $$

อาจลองเขียนโค้ดและใช้คำสั่งของ SymPy ในการแก้สมการเพื่อหารากของสมการได้ดังนี้

s = symbols('s')
roots = solve( Eq(s**4 + 2*s**3 + 2*s**2 + 2*s + 1,0), s )
print(roots)
# => [-1, -I, I]

รากของสมการซึ่งมีดีกรีของพหุนามเท่ากับ 4 หาได้ดังนี้

$$s_{1,2} = -1, \; s_{3,4} = \pm j, \; j := \sqrt{-1} \\ $$ ซึ่งจะเห็นได้ว่า $s_{1,2}$ มีค่าเป็นจำนวนจริงแต่ซ้ำกัน (repeated roots) และ $s_{3,4}$ เป็นจำนวนเชิงซ้อน (เป็นสังยุคของจำนวนเชิงซ้อนคู่กัน)

ดังนั้นรูปแบบของฟังก์ชันที่เป็นผลเฉลยของสมการเชิงอนุพันธ์ที่เป็นเอกพันธ์ เขียนได้ดังนี้

$$y_h(t) = C_1 e^{-t} + C_2 t\, e^{-t} + C_3 cos(t) + C_4 sin(t)$$ โดยที่ $C_1, ..., C_4$ เป็นค่าคงที่ใด ๆ

ผลเฉลยเจาะจงสำหรับสมการเชิงอนุพันธ์ที่ไม่เป็นเอกพันธ์ ก็มีรูปแบบดังนี้

$$y_p(t) = 2 e^{-2t} \\ $$

ดังนั้นเมื่อได้ผลเฉลยทั้งสองกรณีแล้ว ก็นำมาหาผลรวม

$$\begin{align} y(t) &= y_h(t) + y_p(t) \\ &= C_1 e^{-t} + C_2 t\, e^{-t} + C_3 cos(t) + C_4 sin(t) + 2 e^{-2t} \\ \end{align}$$

หากพิจารณาเงื่อนไขสำหรับค่าเริ่มต้น จะได้ผลเฉลยตามรูปแบบของฟังก์ชันดังนี้

$$y(t)= \frac{1}{2} e^{-t} + 6 t\, e^{-t} - \frac{5}{2} cos(t) - \frac{1}{2} sin(t) + 2 e^{-2t} \\ $$

ตัวอย่างการเขียนโค้ด Python เพื่อหาผลเฉลยเชิงตัวเลข มีดังนี้

import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

# ODE: 
#  D4y + 2*D3y + 2*D2y + 2*D1y + y = 10*exp(-2*t);
# Initial conditions:
#  y(0)==0; D1y(0)==1; D2y(0)==-1; D3y(0)==2

t0 = 0.0   # start point for t
tf = 20.0  # end point for t
h  = 0.1   # step size

def f(t,z):
    dz1 = z[1]
    dz2 = z[2]
    dz3 = z[3]
    dz4 = -2*z[3] - 2*z[2] - 2*z[1] - z[0] + 10*np.exp(-2*t)
    return [dz1,dz2,dz3,dz4]

def f_exact(t):
    return (1/2)*np.exp(-t) + (6*t)*np.exp(-t) \
            - (5/2)*np.cos(t) - (1/2)*np.sin(t) + 2*np.exp(-2*t)

ics = [0,1,-1,2]        # vector of initial values
ts = np.arange(t0,tf,h) # array of time steps
# solve the ODE system numerically 
sol = solve_ivp( f, [t0,tf], ics, t_eval=ts, method='RK45' )

plt.clf()
# plot the numerical solution compared with the analytic solution
plt.plot( ts, sol.y[0], 'o', ts, f_exact(ts) )
plt.grid()
plt.xlabel( 't' )
plt.ylabel( 'y(t)' )
plt.legend( ['numeric','analytic'] )
plt.show()

รูป: กราฟของฟังก์ชันที่เป็นผลเฉลย (ใช้ค่า $h=0.1$)

 

---

▷ ตัวอย่างที่ 4#

ถัดไปเป็นการหาผลเฉลยของสมการเชิงอนุพันธ์ แต่จะลองเขียนโค้ด Python โดยใช้วิธีการที่เรียกว่า Runge-Kutta (4th-order)

สมการเชิงอนุพันธ์อันดับสองต่อไปนี้ $$y''(t) + 2y'(t) + y(t) = e^{-t}, \; y(0)=1, y'(0)=0 \\ $$ มีฟังก์ชัน $y(t)$ ที่เป็นผลเฉลยดังนี้ $$y(t) = \frac{1}{2} e^{-t}(t^2 + 2t + 2)$$

 

**การคำนวณด้วยวิธีการของ Runge-Kutta (4th-order) **#

 

วิธีการของ Runge-Kutta (4th-order) หรือ RK4 จะใช้สำหรับการผลเฉลยเชิงตัวเลขของสมการเชิงอนุพันธ์อันดับหนึ่ง และมีการกำหนดค่าเริ่มต้นตามรูปแบบต่อไปนี้ $$y' = f(t,y), \;\;y(0) = y_0 \\ $$

และการคำนวณจะเริ่มจากค่าเริ่มต้น $(t_0,y_0)$ เพื่อคำนวณค่าของ $y$ ในตำแหน่งถัดไปคือ $y_i=y_1, y_2, ...$ สำหรับ $t_i=t_1, t_2, ...$ ตามลำดับ และให้มีระยะห่างคงที่เท่ากับ $h > 0$ (เรียกว่า step-size) หรือ $t_{i+1} = t_i + h$

สูตรในการคำนวณสำหรับ RK4 สำหรับ $i=0,1,2,...$ มีดังนี้

$$\begin{align} k_1 &= f( t_i, y_i ) \\ k_2 &= f( t_i + h/2, y_i + k_1 h/2 ) \\ k_3 &= f( t_i + h/2, y_i + k_2 h/2 ) \\ k_4 &= f( t_i + h, y_i + k_3 h ) \\ y_{i+1} &= y_i + (k_1 + 2k_2 + 2k_3 + k_4)h/6 \\ \end{align}$$ ดังนั้นถ้าจะคำนวณค่า $y_{i+1}$ ในแต่ละรอบ จะต้องทราบค่าของ $y_i, k_1,...,k_4$ และ $h$

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

$$\begin{align} K_1 &= F( t_i, Y_i ) \\ K_2 &= F( t_i + h/2, Y_i + K_1 h/2 ) \\ K_3 &= F( t_i + h/2, Y_i + K_2 h/2 ) \\ K_4 &= F( t_i + h, Y_i + K_3 h ) \\ Y_{i+1} &= Y_i + (K_1 + 2K_2 + 2K_3 + K_4)h/6 \\ \end{align}$$ โดยที่ $Y$ และ $K_1,...,K_4$ เป็นเวกเตอร์ เช่นเดียวกับ $F(t,Y)$ ที่ให้ค่าเป็นเวกเตอร์

ดังนั้นถ้าจะเขียนโค้ด Python สำหรับตัวอย่างที่ 4 โดยใช้ RK4 ก็ทำได้ดังนี้

import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
import math

t0 = 0.0   # start point for t
tf = 10.0  # end point for t
h  = 0.2   # step size

def F(t,Y):
    dY1 = Y[1]
    dY2 =-2*Y[1] - Y[0] + math.exp(-t)
    return np.array([dY1,dY2])

ics = [1.0,0.0]
ts = np.arange(t0,tf,h) # array of time steps

Y = np.zeros( (len(ts),2) ); 
Y[0,:] = np.array(ics)

for i in range(0,len(ts)-1):  
    K1 = F( ts[i]       , Y[i,:]          );
    K2 = F( ts[i] + h/2 , Y[i,:] + h*K1/2 );
    K3 = F( ts[i] + h/2 , Y[i,:] + h*K2/2 );
    K4 = F( ts[i] + h   , Y[i,:] + h*K3   );
    Y[i+1,:] = Y[i,:] + h*(K1 + 2*K2 + 2*K3 + K4)/6; 

f_exact = lambda t: (1/2)*np.exp(-t)*(t**2 + 2*t + 2);

Y = Y.T # transpose
plt.clf()
plt.plot( ts, Y[0], 'o', ts, f_exact(ts), '-' )
plt.grid()
plt.xlabel( 't' )
plt.ylabel( 'y(t)' )
plt.legend(['numeric','analytic'])
plt.show()

รูป: กราฟของฟังก์ชันที่เป็นผลเฉลยจากการใช้วิธี RK4 (ใช้ค่า $h=0.2$)

 

---

▷ กล่าวสรุป#

บทความนี้ได้นำเสนอตัวอย่างการเขียนโค้ดสำหรับ Python เพื่อหาผลเฉลยของสมการเชิงอนุพันธ์สามัญ ด้วยวิธีการคำนวณเชิงสัญลักษณ์ และเปรียบเทียบกับวิธีการคำนวณเชิงตัวเลข

 

---

This work is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License.

Created: 2022-04-14 | Last Updated: 2022-04-29