Добавьте сопротивление воздуха движению снаряда

С силой сопротивления$- \alpha \left|{\dot{\bf r}}\right| {\dot{\bf r}}$и гравитационная сила$-mg {\bf {\hat{y}}}$, уравнения движения (см. мой ответ на этот вопрос )

$$\begin{align} {\ddot{x}} &= - \beta {\dot{x}} \sqrt{{\dot x}^2+{\dot y}^2} \\ {\ddot{y}} &= - g - \beta {\dot{y}} \sqrt{{\dot x}^2+{\dot y}^2} \end{align}$$ где$\beta = \alpha / m$,$\alpha = \rho C_d A / 2$,$\rho$это плотность воздуха,$C_d$- коэффициент лобового сопротивления, а$A$- площадь поперечного сечения снаряда.

Я никогда не видел, чтобы приведенные выше уравнения решались аналитически. Вот некоторый код Python, который отображает решение с сопротивлением воздуха и без него. Для запуска скопируйте в текстовый файл myFile.py и запустите

python myFile.py

из командной строки.

from pylab import *
import math

# Physical constants
g = 9.8
m = 1.0
rho = 1.0
Cd = 1.0
A = math.pi * pow(0.01,2.0)
alpha = rho * Cd * A / 2.0
beta = alpha / m

# Initial conditions
X0 = 1.0
Y0 = 0.0
Vx0 = 70.0
Vy0 = 80.0

# Time steps
steps = 1000
t_HIT = 2.0*Vy0/g
dt = t_HIT / steps

# No drag
X_ND = list()
Y_ND = list()

for i in xrange(steps+1):
  X_ND.append(X0 + Vx0 * dt * i)
  Y_ND.append(Y0 + Vy0 * dt * i - 0.5 * g * pow(dt * i,2.0))

# With drag
X_WD = list()
Y_WD = list()
Vx_WD = list()
Vy_WD = list()

for i in xrange(steps+1):
  X_ND.append(X0 + Vx0 * dt * i)
  Y_ND.append(Y0 + Vy0 * dt * i - 0.5 * g * pow(dt * i,2.0))

# With drag
X_WD = list()
Y_WD = list()
Vx_WD = list()
Vy_WD = list()

X_WD.append(X0)
Y_WD.append(Y0)
Vx_WD.append(Vx0)
Vy_WD.append(Vy0)

stop = 0
for i in range(1,steps+1):
  if stop != 1:
    speed = pow(pow(Vx_WD[i-1],2.0)+pow(Vy_WD[i-1],2.0),0.5)

    # First calculate velocity
    Vx_WD.append(Vx_WD[i-1] * (1.0 - beta * speed * dt))
    Vy_WD.append(Vy_WD[i-1] + ( - g - beta * Vy_WD[i-1] * speed) * dt)

    # Now calculate position
    X_WD.append(X_WD[i-1] + Vx_WD[i-1] * dt)
    Y_WD.append(Y_WD[i-1] + Vy_WD[i-1] * dt)

    # Stop if hits ground
    if Y_WD[i] <= 0.0:
      stop = 1

# Plot results
plot(X_ND, Y_ND)
plot(X_WD, Y_WD)
show()