Модель с рогаткой Python - PullRequest
0 голосов
/ 29 мая 2018

Я пытаюсь смоделировать снаряд из рогатки.

Это мой код:

from pylab import * 
import numpy as np
from scipy.integrate import odeint
import seaborn
## set initial conditions and parameters
g = 9.81            # acceleration due to gravity
th = 30            # set launch angle
th = th * np.pi/180.   # convert launch angle to radians
v0 = 10.0           # set initial speed
c = 0.5             # controls strength of air drag
d = 0.02 # diameter of the spherical rock
A = pi * (d/2)**2 #the cross-sectional area of the spherical rock
ro = 1.2041 #the density of the medium we are perfoming the launch in
m = 0.01 #mass

x0=0                # specify initial conditions
y0=0
vx0 = v0*sin(th)
vy0 = v0*cos(th)

## defining our model
def slingshot_model(state,time):
    z = zeros(4)    # create array to hold z vector
    z[0] = state[2] # z[0] = x component of velocity
    z[1] = state[3] # z[1] = y component of velocity
    z[2] = - c*A*ro/2*m*sqrt(z[0]**2 + z[1]**2)*z[0]         # z[2] = acceleration in x direction
    z[3] = -g/m - c*A*ro/2*m*sqrt(z[0]**2 + z[1]**2)*z[1]       # z[3] = acceleration in y direction
    return z

## set initial state vector and time array
X0 = [x0, y0, vx0, vy0]        # set initial state of the system
t0 = 0
tf = 4 #final time
tau = 0.05 #time step


# create time array starting at t0, ending at tf with a spacing tau
t = arange(t0,tf,tau)   

## solve ODE using odeint
X = odeint(slingshot_model,X0,t) # returns an 2-dimensional array with the 
                        # first index specifying the time and the
                        # second index specifying the component of
                        # the state vector

# putting ':' as an index specifies all of the elements for
# that index so x, y, vx, and vy are arrays at times specified 
# in the time array
x = X[:,0]  
y = X[:,1] 
vx = X[:,2] 
vy = X[:,3]

plt.rcParams['figure.figsize'] = [10, 10]
plot(x,y)

Но это дает мне этот сюжет, который не имеет смысла для меня:

enter image description here

Чего мне не хватает?Ценности не должны выходить, как они, но для жизни я не понимаю, почему.

Это, вероятно, что-то тривиальное, но я слишком долго смотрел на это, поэтому я решил принестив свежем свете глаз лучший курс действий.

1 Ответ

0 голосов
/ 29 мая 2018

Я думаю, что есть как минимум две основные проблемы с вашими вычислениями:

  1. Обычно угол определяется относительно оси X.Поэтому

    vx0 = v0*cos(th) # not sin
    vy0 = v0*sin(th) # not cos
    
  2. Самое главное, почему вы делите ускорение свободного падения g на массу?(см. z[3] = -g/m...) Это не имеет смысла для меня.НЕ делить по массе!

РЕДАКТИРОВАТЬ:

Исходя из вашего комментария и связанных формул, становится ясно, что ваш код также страдает третьей ошибкой: условия воздушного сопротивления должны быть обратно пропорциональны массе :

enter image description here

...