Running Python script, 'Memory error', 'for' loops

Question

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?

You have a double-nested loop to N, so your arrays will be O(N^2) in size. — Tom Dalton
– Tom Dalton, Commented Jan 27, 2016 at 22:33
Thanks, but is there any other way to transform polar to cartesian coordinates for each value of t? I already tried 'pol2cart' but it didn't work. — ffalk
– ffalk, Commented Jan 27, 2016 at 22:39
Is it intentional that plt.plot is inside the for j in ... loop? Don't you want to finish computing x and y before plotting them? — Nick White
– Nick White, Commented Jan 27, 2016 at 22:43
Sorry, N is the dimension of a vector and is set to 100000, while T=300 and dt=T/N. But I've tried running the code with other values — ffalk
– ffalk, Commented Jan 27, 2016 at 22:44

Tom Dalton · Accepted Answer · 2016-01-27 22:57:39Z

1

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!

answered Jan 27, 2016 at 22:57

Tom Dalton

6,19027 silver badges36 bronze badges

Sign up to request clarification or add additional context in comments.

2 Comments

ffalk Over a year ago

I just ran the code using spyder, and while it returned no error, it only plotted one of the odes' solutions. (I've been examining the behaviour manually with various initial conditions, but tried to work out a way to optimize the code in order to make it evaluate different set of IC's for rho[0], and then, depending on those results, vary the range and step of i)

Tom Dalton Over a year ago

It might be worth printing out the first few pairs for x and y, for each i. You might see the values are crazy, or that it's not making it to the second iteration of the i loop, or something else.

Collectives™ on Stack Overflow

Running Python script, 'Memory error', 'for' loops

1 Answer 1

2 Comments

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

2 Comments

Your Answer

Sign up or log in

Post as a guest

Related