使用 Python 建模拋射運動


讓我們首先了解一些基本概念和方程式,這些概念和方程式是拋射運動建模的基礎。下圖解釋了拋射運動的一些基本術語。

這裡,“u”是拋射體被丟擲的速度。𝛼是拋射體的投射角。拋射體在飛行過程中所走的路徑。射程是拋射體水平方向上運動的最大距離。$\mathrm{h_{max}}$是拋射體達到的最大高度。此外,拋射體運動到射程所需的時間稱為飛行時間。

拋射運動分為三類:

  • 水平拋射。

  • 斜拋運動,落在水平地面上。

  • 斜拋運動,但落點在不同高度。

在本教程中,我們將涉及這三種情況。

情況 1:水平拋射

在第一種情況下,拋射體從高處以水平速度(𝑢)發射,並在一段時間後落到地面,落點距離為 R,如下圖所示。

根據運動方程,h、R 和 t(飛行時間)之間的關係可以寫成:

$$\mathrm{h=\frac{1}{2}gt^{2}}$$

$$\mathrm{R=ut}$$

為了用 Python 建模,我們首先需要了解給定的條件和需要求解的結果。為此,我們將從一個問題開始。

示例 1

假設一架飛機在 900 米高空以 140 米/秒的勻速水平速度投下一包食物。這包食物將落在水平距離多遠的地方?

在這個問題中,給出了 ℎ 和 𝑢,需要求解 𝑅。程式如下:

首先,使用 input() 函式從使用者處獲取資料,併為重力加速度定義一個變數。由於鍵盤輸入儲存為字串,因此在 intput() 前面加上 float 進行型別轉換。

# Input data
h=float(input("Enter the height of object\n"))
u=float(input("Enter initial velocity in m/s \n"))
# acceleration due to gravity
g=9.81

根據方程 1 求解 "t",根據方程 2 求解 "R"。

# Time of flight
t=sqrt(2*h/g)
# Range of projectile
R=u*t

但是 sqrt() 函式在常規 Python 包中沒有定義,因此需要在程式開始時呼叫 math 模組。

from math import *

然後列印飛行時間和射程的值作為答案。這裡我們使用了 f-string(格式化字串)來簡化列印過程。

# Output
print(f"Time of flight = {t} s")
print(f"Range = {R} m")

輸出

最後,程式輸出將如下所示:

Enter the height of object
900
Enter initial velocity in m/s
140
Time of flight = 13.545709229571926 s
Range = 1896.3992921400697 m

如果需要,可以使用 round 函式將輸出四捨五入到小數點後幾位,如下所示:

# Output
print(f"Time of flight = {round(t,3)} s")
print(f"Range = {round(R,3)} m")

如果要繪製軌跡,則需要先將射程 (𝑅) 分割成 x 陣列,並獲得拋射體在不同高度 (𝑦) 的位置,然後可以繪製這些資料以獲得球的軌跡。這將在第二種情況下進行說明。

情況 2:斜拋運動,落在水平地面上

在這種情況下,拋射體從水平地面發射,併到達與發射點相同的高度,如下圖所示。

在這種情況下,可以使用以下公式獲得最大高度 $\mathrm{h_{max}}$、射程 (R) 和總飛行時間 (𝑇):

$$\mathrm{h_{max}=\frac{u^{2}\sin^{2}\left ( \alpha \right )}{2g}}$$

$$\mathrm{R=\frac{u^{2}\sin(2\alpha)}{g}}$$

$$\mathrm{T=\frac{2\:u\:\sin\left ( \alpha \right )}{g}}$$

重點是如何求解拋射體軌跡 $\mathrm{y(x)}$。假設在某個時間 "t",拋射體位於點 a,座標為 (𝑥,𝑦),如下圖所示。

然後可以透過以下推導得到粒子位置 "y" 關於 "x" 的函式。

"y" 方向的運動方程:

$$\mathrm{y=u_{y}t+\frac{1}{2}a\:t^{2}=u\sin(\alpha )t-\frac{1}{2}gt^{2}}$$

"x" 方向的運動方程:

$$\mathrm{x=u\cos\left ( \alpha \right )t}$$

$$\mathrm{y=x\tan(\alpha)-\frac{1}{2}\frac{gx^{2}}{u^{2}\cos^{2}\left ( \alpha \right )}}$$

方程 6 是我們可以輕鬆地對拋射運動進行建模並繪製軌跡的方程。

示例 2

為了進行建模,我們需要陣列,因此將使用NumPy模組,此外,由於需要進行繪圖,因此我們將使用pylab,它實際上是matplot.pylab模組。

下面顯示了詳細的程式碼及其說明:

首先呼叫並匯入numpyplyab模組。

# Importing modules
from numpy import *
from pylab import *

然後從使用者處獲取 u 和 $\alpha$。(不要忘記使用 radians() 函式將角度轉換為弧度)

α=float(input("Enter the angle in degrees\n"))
# converting angle in radians
α=radians(α)

# Initial velocity of projectile
u=float(input("Enter initial velocity in m/s \n"))

然後根據公式 (3) 和 (4) 計算拋射體的射程和最大高度。

# Evaluating Range
R=u**2*sin(2*α)/g


# Evaluating max height
h=u**2*(sin(α))**2/(2*g)

現在使用 linspace() 函式建立 "x" 陣列。linspace 接受 3 個引數,即起始點、結束值以及在開始和結束之間建立的點數,並將其考慮在內。

# Creating an array of x with 20 points
x=linspace(0, R, 20)

然後求解公式 6,它將針對不同的 "x" 計算 "y"。

# Solving for y
y=x*tan(α)-(1/2)*(g*x**2)/(u**2*(cos(α))**2 )

然後使用 plot() 函式繪製資料,並進行一些裝飾以使圖形更有意義。如果需要,可以使用savefig()函式儲存圖形。

這裡我們將圖形編號設定為 1,為了提高圖形的可視性,使用函式figure()將 dpi 設定為 300。

# Data plotting
figure(1,dpi=300)
plot(x,y,'r-',linewidth=3)
xlabel('x')
ylabel('y')
ylim(0, h+0.05)
savefig('proj.jpg')
show()

對於 𝑢=5 米/秒和 𝛼=30∘,上述程式的輸出如下所示:

如果需要以表格形式列印資料,則以下程式片段將非常有用。

print(f"{'S. No.':^10}{'x':^10}{'y':^10}")
for i in range(len(x)):
   print(f"{i+1:^10}{round(x[i],3):^10}{round(y[i],3):^10}")

示例

完整的程式程式碼如下所示:

# Importing modules
from numpy import *
from pylab import *

α=float(input("Enter the angle in degrees\n"))
# converting angle in radians
α=radians(α)

# Initial velocity of projectile
u=float(input("Enter initial velocity in m/s \n"))

# Evaluating Range
R=u**2*sin(2*α)/g

# Evaluating max height
h=u**2*(sin(α))**2/(2*g)

# Creating array of x with 50 points
x=linspace(0, R, 20)

# Solving for y
y=x*tan(α)-(1/2)*(g*x**2)/(u**2*(cos(α))**2 )
# Data plotting

figure(1,dpi=300)
plot(x,y,'r-',linewidth=3)
xlabel('x')
ylabel('y')
ylim(0, h+0.05)
savefig('proj.jpg')
show()

print(f"{'S. No.':^10}{'x':^10}{'y':^10}")
for i in range(len(x)):
   print(f"{i+1:^10}{round(x[i],3):^10}{round(y[i],3):^10}")

我們還可以檢查當拋射體以不同角度發射時會發生什麼,然後繪製它們的軌跡。這將非常有助於理解如何將程式碼打包成函式。

在 Python 中,使用 lambda 可以定義單行函式,如下所示:

f= lambda {arguments}: {output statement}

但是,正如您所看到的,這裡有很多行程式碼,因此這將不起作用,而是我們將使用 def 進行編寫,其主體如下:

def (arguments):
   statements1
   statement 2
   .
   .
   .

這裡,我們將從使用者處獲取 u 後面的所有語句打包成一個名為 Projectile_Plot 的函式。該函式接受兩個引數,即 u 和 $\alpha$,如下所示:

# Projectile Function
def Projectile_Plot(u,α):

   # for pronging legend
   αd=α
   
   # converting angle in radians
   α=radians(α)
   
   # Evaluating Range
   R=u**2*sin(2*α)/g
   
   # Evaluating max height
   h=u**2*(sin(α))**2/(2*g)
   
   # Creating array of x with 20 points
   x=linspace(0, R, 20)
   
   # Solving for y
   y=x*tan(α)-(1/2)*(g*x**2)/(u**2*(cos(α))**2 )
   
   # Data plotting
   figure(1,dpi=300)
   plot(x,y,'-',linewidth=2,label=f'α = {αd}')
   xlabel('x')
   ylabel('y')
   ylim(0, h+0.05)
   legend()
   savefig('proj2.jpg')

最後,對於 5 米/秒的速度,可以使用以下程式碼繪製拋射運動(這裡您應該理解,由於角度將發生變化,因此需要使用 for 迴圈):

# Input data
# Initial velocity of projectile
u=float(input("Enter initial velocity in m/s \n"))
for i in range(0,100,10):
Projectile_Plot(u,i)
show()

輸出

程式輸出如下所示:

情況 3:斜拋運動,但落點在不同高度

在本節中,我將向您展示如何在拋射體的起飛點和落地點不在同一高度時對拋射體進行建模,如下圖所示。

在這種情況下,可以看出投射點 (A) 比落地點 (B) 高 $𝑦_{0}$ 個單位。

分析方法如下:

應用垂直方向上的運動方程:

$$\mathrm{s=u_{y}t+\frac{1}{2}a_{y}t^{2}}$$

但是 $\mathrm{u_{y}=u\sin(\alpha)}$ 且 $a_{y}=-g$(即重力加速度)。因此,方程 7 變為

$$\mathrm{y=u\sin(\alpha)t-\frac{1}{2}gt^{2}}$$

到達點 B 所需的時間需要如下計算:

當拋射體到達 B 時,y 將被替換為 $−𝑦_{0}$,因為此時落點低於 X 軸。因此,方程 8 變為:

$$\mathrm{-y_{0}=u\sin(\alpha)t-\frac{1}{2}gt^{2}}$$

為了更好地理解,可以將其簡化為:

$$\mathrm{gt^{2}-2u\sin(\alpha)t-2y_{0}=0}$$

現在需要求解此二次方程(即方程 10)以獲得 t。

最後,水平射程可以計算為:

$$\mathrm{R=u\cos(\alpha)t}$$

在這種情況下,由於需要求解二次方程(即方程 10),因此為了簡化任務,我們將使用 NumPy(數值 Python)。在我們對拋射運動進行程式設計之前,讓我向您展示如何使用 NumPy 求解二次方程的根。

NumPy 中方程的根

假設我們想要找到以下方程的根:

$$ax^{2}+bx+c=0$$

第一步是呼叫 NumPy,如下所示:

from numpy import *

然後將上述多項式的係數儲存在陣列中,如下所示(注意:從多項式的最高次冪到最低次冪):

coeff=array([a,b,c])

然後呼叫函式 roots() 來求解根,如下所示:

r1,r2=roots(coeff)

由於這是一個二次方程,因此我們已經知道有兩個根,因此我們將它們儲存在變數 "r1" 和 "r2" 中。但是,如果多項式的次數增加,則需要相應地增加根變數的數量。下面顯示的示例將闡明上述過程。

求解多項式 x2 - 5x + 6 = 0 的根

# Importing NumPy
from numpy import *

# The polynomial is x^2 -5x+6=0
# The coefficients of the above equation are:
a=1
b=-5
c=6

# Creating array of coefficients
coeff=array([a,b,c])

# Calling function to find the roots
r1,r2=roots(coeff)

# Printing result
print(f"r1 = {round(r1,3)} and r2 = {round(r2,3)}")

執行後,將生成以下輸出

r1 = 3.0 and r2 = 2.0

尋找多項式x3- 2x2- 5x + 6 = 0的根

# Importing NumPy
from numpy import *

# The polynomial is x^3-2x^2-5x+6=0

# The coefficients of the above equation are:

a=1
b=-2
c=-5
d=6

# Creating array of coefficients
coeff=array([a,b,c,d])

# Calling function to find the roots
r1,r2,r3=roots(coeff)

# Printing result
print(f"r1 = {round(r1,3)}, r2 = {round(r2,3)} and r3 = {round(r3,3)}")

執行後,您將獲得以下輸出 -

r1 = -2.0, r2 = 3.0 and r3 = 1.0

現在我們能夠解決斜拋運動問題了。

示例3

一顆子彈從距地面100米高的地方發射,速度為100米/秒,角度為45°。需要評估以下內容 -

  • 飛行時間

  • 飛行距離

  • 拋射體最大高度

解決方案

匯入包 -

from numpy import *
from pylab import *

提供輸入資料 -

# Velocity at the start
u=100
# Angle of projectile at start
α=radians(45)
# Elevation of starting point above the ground
y0=100
# Acceleration due to gravity
g=9.81

計算飛行時間 -

# Time of flight
a=g
b=-2*u*sin(α)
c=-2*y0

# coefficient array
coeff=array([a,b,c])

# finding roots
t1,t2=roots(coeff)
print(f"t1= {t1} and t2= {t2}")

這將返回兩個根,即時間的兩個值。將接受非負值。

拋射體的最大高度將透過以下公式獲得 -

$$\mathrm{h_{max}=\frac{u^{2}\sin^{2}(\alpha)}{2g}+y_{0}}$$

# From the throwing point
h1=u**2*(sin(α))**2/(2*g)

# total
h_max=h1+y0
print(f"h_max= {round(h_max,3)} m")

射程將透過公式11計算,如下所示 -

R=u*cos(α)*max(t1,t2)
# max(t1,t2) will return the positive value

print(f"R= {round(R,3)} m")

因此,完整程式和程式輸出如下所示 -

# Importing packages
from numpy import *
from pylab import *

# Input data
# Velocity at the start
u=100

# Angle of projectile at start
α=radians(45)

# Elevation of starting point above the ground
y0=100

# Acceleration due to gravity
g=9.81

# Time of flight
a=g
b=-2*u*sin(α)
c=-2*y0

# coefficient array
coeff=array([a,b,c])

# finding roots
t1,t2=roots(coeff)
print(f"t1= {t1} and t2= {t2}")

# Maximum Height
# From the throwing point
h1=u**2*(sin(α))**2/(2*g)

# total
h_max=h1+y0
print(f"h_max= {round(h_max,3)} m")

# Range
R=u*cos(α)*max(t1,t2)
#max(t1,t2) will return the positive value

print(f"R= {round(R,3)} m")

執行程式程式碼時,它將產生以下輸出 -

t1= 15.713484026367722 and t2= -1.2974436352046743
h_max= 354.842 m
R= 1111.111 m

繪製拋射運動軌跡

現在任務是繪製此拋射運動軌跡。方程式將與我之前使用的相同,即:

$$\mathrm{y=x\tan(\alpha)-\frac{1}{2}\frac{gx^{2}}{u^{2}\cos^{2}(\alpha)}}$$

繪圖過程如下 -

  • 首先匯入pylab

  • 考慮(0,0)和$(R,-y_{0})$,繪製斜面。

  • 然後從0到R建立x的陣列。

  • 求解公式6以獲得對應於x的y陣列。

  • 最後使用plot()函式繪製(𝑥,𝑦)。

完整的Python程式碼和繪圖如下所示 -

from pylab import *
# Figure name
figure(1, dpi=300)

# Plotting inclined surface
plot([0,R],[0,-y0],linewidth=5)

# plotting y=0 line
plot([0,R],[0,0],'k',linewidth=1)

# Array of x
x=linspace(0,R,50)

# Evaluating y based on x
y=x*tan(α)-(1/2)*(g*x**2)/(u**2*(cos(α))**2 )

# Plotting projectile
plot(x,y,'r-',linewidth=2)
xlabel('x')
ylabel('y')
savefig("Inc_Proj.jpg")
show()

它將產生以下繪圖 -

需要注意的是,只有在評估了拋射體的射程後才能繪製拋射體。因此,您將首先評估射程,然後使用程式碼進行繪製。

更新於: 2023年3月14日

7K+瀏覽量

開啟您的職業生涯

透過完成課程獲得認證

開始學習
廣告