Как решить дифференциальное уравнение 2-го порядка для движения снаряда с сопротивлением воздуха? - PullRequest
Новая
       98

Как решить дифференциальное уравнение 2-го порядка для движения снаряда с сопротивлением воздуха?

0 голосов
/ 21 апреля 2020

уравнение: 

d^2 r/dt^2 = -c/m (dr/dt)+g

, где r - положение снаряда, c - коэффициент сопротивления, m - масса снаряда, g - ускорение под действием силы тяжести.

Предполагая, что только два измерения в форме компонента, это, конечно, гласит: 

d^2 X/dt^2 = -c(dX/dt)= U
d^2 Y/dt^2 = -c/m(dY/dt)+g

, если мы используем методологию, описанную выше, и явно определяем скорости в X И Y как, 

U = dX/dt

и 

V = dX/dt

, тогда вся связанная система уравнений равна, 

dU/dt= -c/m(U)
dV/dt= - c/m(V)+g
dX/dt= U
dY/dt = V

Параметры для этой системы ODE: c = 0,5 кгс ^ - 1, m = 2 кг и g = −9,81 мс ^ -2.

инициализация переменных как (U0, V0, X0, Y0) = (173, 100, 0, 0), которая запускает снаряд из начала координат под углом degrees 30 градусов от горизонтали.

как мне написать новую функцию в python, используя rk4 (я хочу знать, как это закодировать), которая реализует систему из четырех ODE выше, которая решает проблему движения двумерного снаряда ....? Пожалуйста, помогите, я очень плохо знаком с ODE и CODING. СПАСИБО

У меня до сих пор есть следующее ... и оно не работает, и я действительно не знаю, что делать для этой конкретной проблемы c, я должен получить график снаряда а также ... может кто-нибудь, пожалуйста, улучшите мой код, пожалуйста, спасибо 

import numpy as np
import matplotlib.pyplot as plt


def projectileMotion_V(t, M, g, c):
   return -c/M * V0 + g 


def projectileMotion_U(t, c, M):
   return -c/M * U0

V0 = 100         
U0 = 173
ang = 30.0      
c = 0.5       
dt = 0.1  
M = 2.0         
g = -9.81 
h = 0.1

t = [0]                         
x = [0]                         
y = [0]
vx = [V0*np.cos((ang*np.pi)/180)]  
vy = [U0*np.sin((ang*np.pi)/180)]
ax = [-(c*V0*np.cos((ang*np.pi)/180))/M]          
ay = [g-(c*U0*np.sin((ang*np.pi)/180))/M]



def solveODEsWithR4Method(t, x, y, vx, vy, ax, ay):
   t.append(t[0]+dt)                
   vx.append(vx[0]+dt*ax[0])  
   vy.append(vy[0]+dt*ay[0])
   x.append(x[0]+dt*vx[0])    
   y.append(y[0]+dt*vy[0])    
   vel = np.sqrt(vx[0+1]**2 + vy[0+1]**2)   
   drag = c*vel                                    
   ax.append(-(drag*np.cos(ang/180*np.pi))/M)     
   ay.append(-g-(drag*np.sin(ang/180*np.pi)/M)) 
return -c/M * V0 + g 


fig,ax = plt.subplots()
ax.plot(t, M, g, c)
plt.show()

Ответы [ 2 ]

1 голос
/ 21 апреля 2020

Как больше ничего не происходит в эти времена пандемии c заблуждений о phantasmi c вирусах-убийцах ...

Пожалуйста, не называйте это сопротивление воздуха, то есть конкретно k*|v|*v. Поскольку для того, чтобы это была сила, коэффициент k должен был бы иметь единицы kg/m, а это не то, что вам дано, ваши формулы для сопротивления, вероятно, верны. Вместо этого назовите это «среднее сопротивление», сопротивление воды будет таким. векторные состояния 

def RK4step(f,u,dt):
    k1 = dt*f(u)
    k2 = dt*f(u+0.5*k1)
    k3 = dt*f(u+0.5*k2)
    k4 = dt*f(u+k3)
    return u + (k1+2*k2+2*k3+k4)/6

def RK4integrate(f, u0, tspan):
    u = np.zeros([len(tspan),len(u0)])
    u[0,:]=u0
    for k in range(1, len(tspan)):
        u[k,:] = RK4step(f, u[k-1], tspan[k]-tspan[k-1])
    return u

и применение обоих кодов вместе для вычисления траектории 

dt = .1
t = np.arange(0,10,dt)
u0 = np.array([0, 0, 173, 100])

sol_RK4 = RK4integrate(motion, u0, t)
x,y,vx,vy = sol_RK4.T
plt.plot(x,y)
0 голосов
/ 24 апреля 2020

Set v = dr/dt. Тогда уравнения становятся: 

dv/dt = - (c/m) * v  +  g

Умножим обе части уравнения на exp((c/m)*t): 

exp((c/m)*t) * dv/dt = - exp((c/m)*t) * (c/m) * v  +  exp((c/m)*t) * g

Переместим первый столбец слагаемых с правой стороны на левую : 

exp((c/m)*t) * dv/dt + exp((c/m)*t) * (c/m) * v  =  exp((c/m)*t) * g

Затем по правилу дифференцирования произведения, примененному к exp((c/m)*t) * v, можно получить 

d/dt( exp((c/m)*t) * v )  =  exp((c/m)*t) * g

d/dt( exp((c/m)*t) * v )  =  d/dt( (m/c) * exp((c/m)*t) * g )

Теперь вы можете интегрировать обе стороны относительно t, и существует вектор u0 такой, что 

exp((c/m)*t) * v  =  u0  +  (m/c) * exp((c/m)*t) * g

Умножим обе стороны на exp(-(c/m)*t): 

v  =  exp(-(c/m)*t) * u0  +  (m/c) * g

Если начальная скорость равна v0, то вам нужно установить u0 = v0 - (m/c)*g, чтобы Окончательная формула для скорости: 

v  =  exp(-(c/m)*t) * ( v0 - (m/c)*g )  +  (m/c) * g

Интегрируйте v еще раз относительно t, и вы получите формулу для позиции: 

r  = r0 - (m/c) * ( v0 + (m/c)*g )  - (m/c) * exp(-(c/m)*t) * ( v0 - (m/c)*g )  +  (m/c) * t * g

, где r0 - начальная позиция.

Итак, окончательные формулы для решения: 

r  = r0 - v0 + (m/c)*g  +  exp(-(c/m)*t) * ( v0 - (m/c)*g )  +  (m/c) * t * g

v  =  exp(-(c/m)*t) * ( v0 - (m/c)*g )  +  (m/c) * g
...