如何求解具有空气阻力的弹丸运动的二阶微分方程?

Question

等式为：

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

该常微分方程系统的参数为 c = 0.5 kgs^−1， m = 2kg 和 g = −9.81 ms^−2。

将变量初始化为(U0, V0, X0, Y0) = (173, 100, 0, 0)，从原点以与水平面 ∼ 30 度的角度发射弹丸。

如何使用 rk4 在 python 中编写一个新函数(我想知道如何编码)，该函数实现了上述四个常微分方程的系统，解决了 2D 弹丸运动问题....？ 请帮忙 我对常微分方程和编码很陌生。谢谢

---

到目前为止，我有以下内容...它不起作用，我真的不知道该怎么做这个特定问题，我也打算获得弹丸图......有人可以改进我的代码吗，谢谢

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

Answer 1

由于在这个关于幻觉杀手病毒的妄想大流行的时代没有发生任何其他事情......

请不要称其为空气阻力，这是特别k*|v|*v。至于这是一个力，系数k需要有单位kg/m，这不是你得到的，你的电阻公式可能是正确的。称其为"中等阻力"，防水性就是这样。

然后对加速进行编码

c = 0.5; m = 2; g = -9.81;
def motion(x,v):
x,y,vx,vy = v
return np.array([vx,vy, -c/m * vx + g, -c/m * vy ])

从面向向量状态的地方复制 RK4 代码

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)

Answer 2

设置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