Запуск сценария Python, «Ошибка памяти», «для» циклов

Question

Я пытаюсь смоделировать проблему с тремя телами, где у меня есть три обобщенных координаты (один радиальный и два угловых) и три дифференциальных уравнения второго порядка (связанные) , Я хочу видеть, как система развивается при изменении начальных условий, rho [0]. Что-то не так в моем скрипте? Переменные и параметры хорошо определены выше, поэтому я их опускаю; вот код:

for i in range(1000000,10000000,1000000): 
     rho[0] = i 

     rho[1] = rho[0] + vrho*dt 
     theta[1] = theta[0] + vtheta*dt #angulos radianes 
     phi[1] = phi[0] + vphi*dt 

     for t in range(1, N-1): 
     #Velocidades de las coordenadas 
      v[0,t-1] = (rho[t] - rho[t-1])/dt 
      v[1,t-1] = (theta[t] - theta[t-1])/dt 
      v[2,t-1] = (phi[t] - phi[t-1])/dt 

     #"Ecuaciones diferenciales" 
      rho[t+1] = (2*rho[t] - rho[t-1]) + (rho[t]*(v[2,t-1]**2) - G*M/(rho[t]**2) - (9*G*M*((np.cos(theta[t]-phi[t]))**2)*(l**2))/(8*(rho[t]**4)))*(dt**2) 
      theta[t+1] = (2*theta[t] - theta[t-1]) + ((-3*G*M*np.sin(2*(theta[t]-phi[t])))/(2*(rho[t]**3)))*(dt**2) 
      phi[t+1] = (2*phi[t] - phi[t-1]) + ((3*G*M*np.sin(2*(theta[t]-phi[t]))*(l**2))/(8*(rho[t]**5)) - (2*v[0,t-1]*v[2,t-1])/(rho[t]))*(dt**2) 

      for j in range(0,N): 
       x[j] = rho[j]*np.cos(phi[j]) 
       y[j] = rho[j]*np.sin(phi[j]) 

       plt.plot(x,y) 

И это ошибка:

--------------------------------------------------------------------------- 
    MemoryError        Traceback (most recent call last) 
<ipython-input-3-4c789e1dfc65> in <module>() 
    21    y[j] = rho[j]*np.sin(phi[j]) 
    22 
---> 23    plt.plot(x,y) 

То, что я пытался сделать, чтобы сделать питон решить ОДУ при различных начальных условиях, сэкономив ро, тета, фи и V вектор в ордене, чтобы манипулировать ими. Учитывая тот факт, что я хочу траекторию тела, для тех значений rho и phi, я хочу превратить их в декартовую, построить траекторию, а затем перезапустить ODE со следующим начальным условием для rho.

Мне пришлось ввести счетчик j, так как x и y были определены как np.zeros (N), чтобы соответствовать размерам rho, theta и phi. Счетчик представляет положение такого вектора, и идея состоит в том, что для каждой позиции в rho и phi декартовый эквивалент сохраняется в том же положении. Так как t поднимается до N-1, и, следовательно, последний t доступен N-2, я не знал, как избежать добавления другого; чтобы решить разницу в исходных точках, я мог бы написать последнее переключение j для t-1, но разве это не приведет к размерной ошибке?

Answer 1

Вы можете попробовать следующий код:

for i in range(1000000,10000000,1000000): 
    rho[0] = i 

    rho[1] = rho[0] + vrho*dt 
    theta[1] = theta[0] + vtheta*dt #angulos radianes 
    phi[1] = phi[0] + vphi*dt 

    for t in range(1, N-1): 
     # Velocidades de las coordenadas 
     v[0,t-1] = (rho[t] - rho[t-1])/dt 
     v[1,t-1] = (theta[t] - theta[t-1])/dt 
     v[2,t-1] = (phi[t] - phi[t-1])/dt 

     # "Ecuaciones diferenciales" 
     rho[t+1] = (2*rho[t] - rho[t-1]) + (rho[t]*(v[2,t-1]**2) - G*M/(rho[t]**2) - (9*G*M*((np.cos(theta[t]-phi[t]))**2)*(l**2))/(8*(rho[t]**4)))*(dt**2) 
     theta[t+1] = (2*theta[t] - theta[t-1]) + ((-3*G*M*np.sin(2*(theta[t]-phi[t])))/(2*(rho[t]**3)))*(dt**2) 
     phi[t+1] = (2*phi[t] - phi[t-1]) + ((3*G*M*np.sin(2*(theta[t]-phi[t]))*(l**2))/(8*(rho[t]**5)) - (2*v[0,t-1]*v[2,t-1])/(rho[t]))*(dt**2) 

    for j in range(0,N): 
     x[j] = rho[j]*np.cos(phi[j]) 
     y[j] = rho[j]*np.sin(phi[j]) 

    plt.plot(x,y) 

Я думаю, что из-за плохой отступа на этапе черчения, то есть каждый раз, когда вы хотели, чтобы построить одну строку, вы на самом деле построения N^2 линии!