使用 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模块。

下面显示了详细的代码及其说明：

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

# 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" 中。但是，如果多项式的次数增加，则需要相应地增加根变量的数量。下面显示的示例将阐明上述过程。

求解多项式 x² - 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³- 2x²- 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+浏览量

开启您的职业生涯

通过完成课程获得认证

开始学习