使用Python建模抛射运动

首先，让我们了解一些基本概念和方程，基于这些概念和方程可以对抛射运动进行建模。下面显示的图解释了抛射运动的一些基本术语。

在这里，”u”是投射物体的速度。𝛼是投射物体的投射角度。弹道是投射物体在飞行过程中所经过的路径。射程是投射物体所行进的最大水平距离。\mathrm{h_{max}}是投射物体所达到的最大高度。此外，投射物体到达射程所花费的时间称为飞行时间。

投射运动可以分为三类−

• 水平投射。

• 倾斜投射在水平地面上。

• 倾斜投射但击中不同高度。

在本教程中，我们将涉及这三种情况。

案例1：水平投射

在第一种情况下，投射物体从一个高处仅具有水平速度(𝑢)，在一定距离R后将击中地面，如下图所示。

根据运动方程，h、R和t（飞行时间）之间的关系可以写成：

h = 1/2 * g * t^2

R = ut

要在Python中模拟它，首先我们必须明确给定和所需的内容。为此，我们将从一个问题开始。

示例1

假设一个飞机在高度900米处以均匀的水平速度140米/秒释放一包食物。该包裹落地的水平距离是多少。

在这个问题中，给定h和u，而需要求解R。程序如下：

首先，使用input()函数从用户那里获取数据，并为重力加速度定义一个变量。由于键盘输入存储为字符串，因此将其强制转换为浮点数，将float放在input()前面。

# 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")

如果想要绘制弹道轨迹，首先需要将射程（R）分成一组x数组，并获得不同高度上的弹丸位置（y），然后可以将这些数据绘制成轨迹。这将在第二种情况下进行解释。

案例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模块。

下面展示了详细的代码及其解释 –

首先调用并导入numpy和pylab模块。

# Importing modules
from numpy import *
from pylab import *

然后从用户那里询问你和\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
   .
   .
   .

在这里，我们将把要求用户输入之后的所有语句捆绑到一个名为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 m/s的抛体运动，可以使用以下代码绘制（在这里您应该理解，由于角度将会变化，因此需要使用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（Numerical 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

寻找多项式x的根 3 - 2x 2 - 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 m/s和45^{0}的速度被发射到离地面100 m的地方。需要评估以下内容：

• 飞行时间

• 飞行距离

• 抛物线的最大高度

解答

引入包 –

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()

它将生成以下图表 −

这里需要注意的是，只有在评估了弹射物的射程之后才能绘制出来。因此，你将首先评估射程，然后使用绘制代码。