Python中带有drag的物理项目符号立即位于z以下

Question

我正在编写一个程序，用不同的参数来计算一颗子弹以20度角射出的距离，不管有没有阻力。出于某种原因，我的程序立即计算出子弹的y值低于零，我想这可能是我使用的物理方程有问题。我没有最扎实的物理学背景，所以在尽我最大的努力找出我哪里出错后，我想我会寻求帮助。在下面的图表中，我试图显示与不同时间步相关的不同距离。你知道吗

    import math
    import matplotlib.pyplot as plt
    rho = 1.2754
    Af = math.pi*(16*0.0254)**2/4
    g = -9.8 #m/s^2
    V = 250.0 #m/s
    #equation of motion
    def update(r,V,a,dt):
        return r+V*dt +.5*a*dt*dt, V + a*dt
    #drag force from velocity vector return x and y 
    def drag(Vx,Vy):
        Vmag = math.sqrt(Vx**2+Vy**2)
        Ff = .5*Cd*rho*Af*Vmag**2
        return -Ff * Vx/Vmag, -Ff *Vy/Vmag

    def getTrajectory(dt, th):
        #initialize problem [time,x,y,Vx,Vy]
        state = [[0.0,0.0,0.0, V*math.cos(th), V*math.sin(th)]]

        while state[-1][2] >= 0 : # while y > 0 
            time = state[-1][0] + dt
            Fd = drag(state[-1][3],state[-1][4]) # Vx,Vy
            ax = Fd[0]/m # new x acceleration
            ay = (Fd[1]/m) +g # new y acceleration

            nextX, nextVx = update(state[-1][1], state[-1][3],ax,dt)                                         #x,Vx,acceleration x, dt
            nextY, nextVy = update(state[-1][2], state[-1][4],ay,dt) #y,Vy, acceleration y, dt
            state.append([time,nextX,nextY,nextVx,nextVy])
            # linear interpolation 
            dtf = -state[-1][2]*(state[-2][0]-state[-1][0])/(state[-2][2]-state[-1][2])

            xf,vxf = update(state[-1][1],state[-1][3], ax, dtf)
            print(xf,vxf)
            yf,vyf = update(state[-1][2],state[-1][4], ay, dtf)
            print(yf,vyf)
            state[-1] = [time+dtf,xf,yf,vxf,vyf]
            return state

    delt = 0.01 #timestep
    th = 20.0 * math.pi/180.0 #changing to radians
    m = 2400 #kg

    Cd = 1.0
    tsteps = [0.001,0.01,0.1,1.0,3.0] # trying different timesteps to find the optimal one
    xdistance = [getTrajectory(tstep,th)[-1][1] for tstep in tsteps]
    plt.plot(tsteps,xdistance)
    plt.show()

一个timestep之后包含time、x、y、Vx和Vy的所有值的状态变量的打印产生： [[0.0, 0.0, 0.0, 234.9231551964771, 85.50503583141717] ，[0.0，0.0，-9.540979117872439e-18234.9231551964771，85.50503583141717]] 第一个数组是初始条件，Vx和Vy被分解成20度角的向量。第二个数组应该让子弹以正的x和y飞走，但是由于某些原因，y值被计算为-9.54e-18，这就放弃了getTrajectory方法，因为它只应该在y>；0时继续

为了简化我的问题，我认为具体研究一下我的update函数可能会有所帮助。你知道吗

    def update(initial,V,a,dt): #current x or y, Velocity in the x or y         direction, acceleration in x or y, timestep
        return initial+V*dt +.5*a*dt*dt, V + a*dt

这个函数应该根据传递的内容返回一个更新的x，Vx或更新的y，Vy。我想这可能是物理方程的问题所在，但我不完全确定。你知道吗

Answer 1

我想我看到了发生了什么：在while循环中，首先向前一步，从它们在t = 0的值计算运动变量，然后把dtf设为负值，然后向后一步，从刚刚在t = 0.001计算的值计算新的t = 0值。我不知道为什么要这样做，但无论如何，你的update()函数使用的是Euler method，这在数值上是不稳定的——当你采取越来越多的步骤时，它往往以某种不可预测的方式越来越偏离真正的解。在本例中，它显示为负y位置，而不是像您所期望的那样将其设置回0。你知道吗

数值积分的任何方法都不会给出完美的结果，所以你总是会期望有一点不精确，可能是正的，也可能是负的。仅仅因为某件事有点消极，就不必担心。也就是说，你可以通过使用一种更稳健的积分方法，比如Runge-Kutta四阶（又称RK4）或辛积分，得到更好的结果。但这些都是一些先进的技术，所以我现在不会太担心它们。你知道吗

至于手头的问题，我想知道你的方程中是否有符号错误，可能是在设置dtf的行中。不过，我不能肯定，因为我真的不明白getTrajectory()应该如何工作。你知道吗

Answer 2

这似乎是一个缩进错误。从线性插值开始确定y过零位置的部分肯定不是积分循环的一部分，尤其不是return语句。修复此缩进级别将提供输出

(3745.3646817336303, 205.84399509886356)
(6.278677479133594e-07, -81.7511764639779)

(3745.2091465034314, 205.8432621147464)
(8.637877372419503e-05, -81.75094436152747)

(3743.6311931495325, 205.83605594288176)
(0.01056407140434586, -81.74756057157798)

(3727.977248738945, 205.75827166721402)
(0.1316764390634169, -81.72155620161116)

(3669.3437231150106, 205.73704583954634)
(10.273469073526238, -80.52193274789204)