我如何在PYTHON中求解带有空气阻力的弹丸运动的二阶微分方程？

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

此ODE系统的参数为c = 0.5 kgs ^ −1，m = 2kg，g = −9.81 ms ^ −2。

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

如何使用rk4在python中编写新函数（我想知道如何对此进行编码），该函数实现了上面包含四个ODE的系统，解决了2D弹丸运动问题...。？请帮助我，我对ODE和编码非常陌生。谢谢

Answer 1

到目前为止，我已经掌握了以下内容……而且它无法正常工作，我真的不知道该如何解决该特定问题，我也打算获得一个弹丸图……有人可以改善我的代码吗？请谢谢

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

'''