Reputation: 59
I'm trying to model the three-body problem, where I have three generalized coordinates (one radial and two angular) and three second order (coupled) differential equations. I want to see how the system evolves while changing the initial conditions, rho[0]. Is there anything wrong in my script? The variables and parameters are well defined above, so I will omit them; here's the code:
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)
And this is the error:
---------------------------------------------------------------------------
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)
What I tried to do was to make python solve the ODE's for different initial conditions, saving rho, theta, phi and the V vector in orden to manipulate them. Given the fact that I want the trajectory of the body, for those values of rho and phi, I want to turn them to cartesian, plot the trajectory, and then restart the ODE's with the next initial condition for rho.
I had to introduce the j counter since x and y were defined as np.zeros(N) in order to match the dimensions of rho, theta and phi. The counter represents the position of such vector, and the idea is that for each position in rho and phi, the cartesian equivalent is saved in the same position. Since t goes up to N-1, and therefore the last t available is N-2, I didn't know how to avoid adding another for; in order to solve the starting point difference, I could write the last 'for' switching j for t-1, but wouldn't that cause a dimensional error?
Upvotes: 2
Views: 878
Reputation: 6190
Can you try the following code:
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)
I think there was some bad indentation at the plotting stage, meaning every time you wanted to plot a single line, you were actually plotting N ^ 2 lines!
Upvotes: 1