Python 2.7 Runge Kutta Orbit GUI

Question

Я пытаюсь создать графический интерфейс в Python, используя Tkinter, который отображает орбитальный путь массы с использованием метода Рунге-Кутты. Мой графический интерфейс работает отлично, но моя проблема заключается в том, что он отображает только прямую линию независимо от того, какие значения я вводил в графический интерфейс. Я надеялся, что кто-то может показать мне, что не так с моей функцией в графическом интерфейсе, так что он будет правильно строить орбиту.

def calcPath(self): 
    M = float(self.entM.get()) 
    m = float(self.entm.get()) 
    G = 6.67e-11 
    c = 3e8 



    velocity = np.array([float(self.entvx.get()),float(self.entvy.get()),float(self.entvz.get())]) 
    pos = np.array([float(self.entx.get()), float(self.enty.get()), float(self.entz.get())]) 

    Force = lambda pos: (G*m*M/(pos**2)) 

    #assigning variables empty lists to append x, y and z values to 
    a = [] 
    b = [] 
    c = [] 
    t = 0.0 
    tf = float(self.enttf.get()) 
    dt = float(self.entdt.get()) 

    while t < tf: # creating a while loop to trace the path for a set period of time 
     #using the Runge-kutta formula from the problem sheet to assign variables at different steps and half steps 
     k1v=(dt*Force(pos))/m 
     k2v=(dt*Force(pos+(k1v/2.0)))/m 
     k3v=(dt*Force(pos+(k2v/2.0)))/m 
     k4v=(dt*Force(pos+k3v))/m 
     velocity=velocity+(k1v/6.0)+(k2v/3.0)+(k3v/3.0)+(k4v/6.0) #calaculating the weighted average of the k values 
     pos=pos+velocity*dt #velocity is not a function of space or time so it will be identical at all 4 k values 

     a.append(pos[0]) # appending the lists for each vaiable to produce a plot 
     b.append(pos[1]) 
     c.append(pos[2]) 
     t = t+dt # setting the time steps 

    #generating the path plot figure and formatting it  
    ax = Axes3D(self.fig) #creating a 3D figure 
    ax.set_title("Path of planetary mass") #produces title on the figure 
    ax.plot3D(a,b,c, color='darkviolet', label='Runge-kutta path method') 
    ax.set_xlabel('x') 
    ax.set_ylabel('y') 
    ax.set_zlabel('z') 
    ax.legend(loc='lower left') #selecting the location of legend 
    self.canvas.show() 

Answer 1

Ваша физика ошибочна. И ваше отношение к уравнению движения неверно.

Сила вычисляет, как

Force = lambda pos: -G*m*M/(norm2(pos)**3)*pos 

и уравнение движения является

м · поз '' (T) = Force (поз (т))

, где

pos '' (t) = d²/dt² pos (t)

, которые вы должны преобразовать в системе первого порядка

[ pos' , vel' ] = [ vel, Force(pos)/m ] 

К этой системе первого порядка теперь вы можете применить метод RK4.

k1p = dt *  vel 
    k1v = dt * Force(pos  )/m 
    k2p = dt *  (vel+0.5*k1v) 
    k2v = dt * Force(pos+0.5*k1p)/m 
    k3p = dt *  (vel+0.5*k2v) 
    k3v = dt * Force(pos+0.5*k2p)/m 
    k4p = dt *  (vel+ k3v) 
    k4v = dt * Force(pos+ k3p)/m 
    pos = pos + (k1p+2*k2p+2*k3p+k4p)/6.0 
    vel = vel + (k1v+2*k2v+2*k3v+k4v)/6.0 

Обратите внимание, как k векторов чередовать между положением и скоростью.

---

Или вы можете сделать корректную реализацию RK4 автоматический с помощью

System = lambda pos, vel : vel, Force(pos)/m 

, а затем внутри цикла времени

k1p, k1v = System(pos   , vel   ) 
    k2p, k2v = System(pos+0.5*k1p*dt, vel+0.5*k1v*dt) 
    k3p, k3v = System(pos+0.5*k2p*dt, vel+0.5*k2v*dt) 
    k3p, k3v = System(pos+ k3p*dt, vel+ k3v*dt) 
    pos = pos + (k1p+2*k2p+2*k3p+k4p)*dt/6.0 
    vel = vel + (k1v+2*k2v+2*k3v+k4v)*dt/6.0 

---

И если вы просто хотите, численные решения, т.е. , не нужно доказывать, что вы можете реализовать RK4, затем используйте scipy.integrate.ode или scipy.integrate.odeint.