จำลองการเคลื่อนที่แบบ Projectile ด้วย Matplotlib [ตอน 1]

[ตอนที่ 2]
การเคลื่อนที่แบบ projectile [1] เป็นรูปแบบการเคลื่อนที่ที่พบได้บ่อยในชีวิตประจำวันโดยเฉพาะในช่วงสงกรานต์ เมื่อมีการฉีดน้ำหรือสาดน้ำก็จะได้เห็นการเคลื่อนที่แบบ projectile เมื่อนั้น

การเคลื่อนที่แบบ projectile เกิดขึ้นเมื่อวัตถุถูกทำให้เคลื่อนที่ในตำแหน่งใกล้กับผิวโลก วัตถุจะเคลื่อนที่เป็นเส้นโค้งพาราโบล่า [4] ก่อนตกลงสู่พื้นด้วยแรงโน้มถ่วง

บทความนี้จะแสดงการใช้ Matplotlib และ Python ในสร้างแบบจำลองการเคลื่อนที่แบบ projectile แต่ก่อนที่จะเข้าสู่การเขียนชุดคำสั่ง ผู้ที่ยังไม่ทราบที่มาที่ไปควรจะทำการศึกษาสมการการเคลื่อนที่ได้จาก [1][2] [3][5] และการใช้ animation module ของ matplotlib [6] เสียก่อน

กิจกรรมที่ต้องคำนวณ
จาก [3] ทำให้เราทราบว่าควรมีชุดคำสั่งที่เกี่ยวข้องกับการเคลื่อนที่ดังนี้

1. การหาค่าการขจัด (displacement) ในแนวตั้งและแนวนอน ของเวลา t ใดๆ
2. การหาค่าระยะสูงสุดในแนวตั้ง (height)
3. หาค่าระยะไกลสุดในแนวนอน (range)
4. หาเวลาในการเดินทางทั้งหมด

ค่าคงที่ (constant variables) ต่าง ๆ ดังนี้
1. gravity (g) มีค่าประมาณ 9.81 m/s²
2. ความเร็วเริ่มต้นของวัตถุ v₀ m/s
3. มุมที่แนวการเคลื่อนที่ของวัตถุต่อพื้นราบ (theta) หน่วยเป็น radian

แปลงความต้องการเหล่านี้ออกมาเป็น Python script


from matplotlib import pyplot as plt
from matplotlib import animation
import numpy as np

เนื่องจากมีการคำนวณค่อนข้างมาก ดังนั้นเราจะใช้ numpy module [7] ซึ่งเป็น module หลักที่นิยมใช้ในการคำนวณในภาษา Python


g = 9.81 #value of gravity
v0 = 10  #initial velocity 
theta = np.radians(45) #angle of launch in radian

กำหนดค่าคงที่ที่จะใช้งานโปรแกรม


def compute_position(v,theta,t):
        sx = v * np.cos(theta) * t
        sy = v * np.sin(theta) * t - 0.5 * g * (t**2)
        return (sx,sy)

ฟังก์ชั่น compute_position รับค่า v (velocity),theta (angle),t (time) เข้ามาเพื่อคำนวณค่าขจัดแนวนอน (sx) และค่าขจัดแนวตั้ง (sy) [3]


def compute_total_time(velocity,theta,gravity):
        t = (2 * velocity * np.sin(theta)) / gravity 
        return t

ฟังก์ชั่น compute_total_time รับค่า velocity (ความเร็ว), theta (angle of launch), gravity (ความเร่ง) เข้าเพื่อคำนวณหาเวลาที่วัตถุจะใช้ในการเดินทางทั้งหมด [3]


def compute_range(velocity,theta,gravity):
        r = (velocity**2) *  np.sin(2 * theta) / gravity
        return r

ฟังก์ชั่น compute_range รับค่า velocity (ความเร็ว), theta (angle of launch), gravity (ความเร่ง) เข้าเพื่อคำนวณหาระยะทางมากที่สุดในแนวนอนในการเคลื่อนที่ [3]


def compute_height(velocity,theta,gravity):
        h = (velocity**2) *  np.sin(2 * theta) / gravity
        return h

ฟังก์ชั่น compute_height รับค่า velocity (ความเร็ว), theta (angle of launch), gravity (ความเร่ง) เข้าเพื่อคำนวณหาระยะสูงสุดในแนวตั้งของการเคลื่อนที่ [3]


t_max = compute_total_time(velocity=v0,theta=theta, gravity=g)
t_index = np.arange(0,t_max,0.01)

x_max = compute_range(velocity=v0,theta=theta,gravity=g)
y_max = compute_height(velocity=v0,theta=theta,gravity=g)

นำมาคำนวณหาพารามิเตอร์ที่ต้องการ ขอให้สังเกตุว่ามีการสร้างตัวแปรชื่อ t_index ขึ้นมาซึ่งเป็น array จะมีค่าตั้งแต่ 0 (เริ่มเคลื่อนที่)จนถึง t_max (วัตถุตกถึงพื้น) ค่าของแต่ละช่วงห่างกัน 0.01 หน่วย (อาจตีความว่าเป็นหน่วยวินาทีก็ได้) ค่าแต่ละค่าใน t_index จะถูกนำส่งไปใช้ในการคำนวณ sx, sy ของแต่ละเฟรม


fig = plt.figure()
ax = plt.axes(xlim=(0, x_max+1), ylim=(0, y_max+1))
ax.set_aspect('equal')

circle1 = plt.Circle((0,0), 0.2, color='red') # create particle
ax.add_patch(circle1)

สร้างพื้นที่ในการสร้างภาพ และสร้าง particle (circle1)


def animate(i,circle): 
        t = t_index[i]
        (sx,sy) = compute_position(v0,theta,t)
        circle.center = (sx,sy)
 
        return circle,

ฟังก์ชั่น animate รับ parameter มาสองตัวคือ i (ลำดับของ frame) และ circle (วัตถุที่ใช้แสดงการเคลื่อนที่) ค่า i จะถูกนำไปหาค่าของ time ใน t_index ที่สอดคล้องกับตำแหน่งของเฟรม แล้วนำไปคำนวณหาตำแหน่งที่สอดคล้องกันใน t_index ค่าที่ได้คือตำแหน่งของวัตถุ (sx,sy) จากนั้นนำค่าที่ได้นี้ไปเปลี่ยนจุดศูนย์กลาง circle ใหม่



anim = animation.FuncAnimation(
 fig=fig, 
 fargs=(circle1,),
 func = animate, 
 frames=len(t_index), 
 interval=30,
 blit=True)

plt.show()

ชุดคำสั่งนี้ทำหน้าที่เช่นเดียวกับที่กล่าวไว้ใน [6]

Complete Script


from matplotlib import pyplot as plt
from matplotlib import animation
import numpy as np

def compute_position(v,theta,t):
 sx = v * np.cos(theta) * t
 sy = v * np.sin(theta) * t - 0.5 * g * (t**2)
 return (sx,sy)

def compute_total_time(velocity,theta,gravity):
 t = (2 * velocity * np.sin(theta)) / gravity 
 return t

def compute_range(velocity,theta,gravity):
 r = (velocity**2) *  np.sin(2 * theta) / gravity
 return r

def compute_height(velocity,theta,gravity):
 h = (velocity**2) * (np.sin(theta)**2)/(2*gravity)
 return h

g = 9.81 #value of gravity
v0 = 10  #initial velocity 
theta = np.radians(45) #angle of launch in degrees

t_max = compute_total_time(velocity=v0,theta=theta, gravity=g)
t_index = np.arange(0,t_max,0.01)


x_max = compute_range(velocity=v0,theta=theta,gravity=g)
y_max = compute_height(velocity=v0,theta=theta,gravity=g)

fig = plt.figure()
ax = plt.axes(xlim=(0, x_max+1), ylim=(0, y_max+1))
ax.set_aspect('equal')

circle1 = plt.Circle((0,0), 0.2, color='red') # create particle
ax.add_patch(circle1)

def animate(i,circle): 
 t = t_index[i]
 (sx,sy) = compute_position(v0,theta,t)
 circle.center = (sx,sy)
 
 return circle,

anim = animation.FuncAnimation(
 fig=fig, 
 fargs=(circle1,),
 func = animate, 
 frames=len(t_index), 
 interval=30,
 blit=True)

plt.show()

ผลลัพธ์ที่ได้