-2

I'm aiming to compute planetary motion using Verlet algorithms.

The problem is, the code takes roughly 30 minutes to do a full cycle. This is basic computation so I have no idea why it's taking so long.

I've read around on here a bit and using numpy should be faster right?

How do I implement it though? Also, it should be coming full circle and stopping at 0 ( or 2*pi*r_earth, if it's calculating distance travelled ), but it doesn't stop at all, just carries on, can anyone see why my limit isn't working?

Here's my code:

grav    =   6.673481*(10**-11)             # = 6.673481e-11  # a direct way
m_sun   =   1.989*(10**30)                 # = 1.989e+30
m_earth =   5.972*(10**24)                 # = 5.972e+24    
r_earth = 149.59787*(10**9)                # = 1.4959787e+11

def verlet_v( prev_v, prev_f, mass ):

    current_v = ( prev_v + ( prev_f / mass ) )
    print "New Velocity: %f" % current_v
    return current_v

def verlet_x( prev_x, current_v ):

    current_x = prev_x + current_v
    print "New Position: %f" % current_x
    return current_x

verlet_v( 20, 50, 3 )

v = 29.8*(10**3)                          # = 2.98e+04
x = 0
f = ( -grav * ( ( m_earth * m_sun ) / r_earth**2 ) )
v_history = []
x_history = []

while( abs(x) > ( -0.1 ) ):               # ABS( <_anything_> ) is ALWAYS > -0.1

    print "Mod(x): %f" % abs(x)
    print "Limit: %f"  % (0)
    v = verlet_v( v, f, m_earth )
    x = verlet_x( x, v )
    v_history.append(v)
    x_history.append(x)

print v_history
print x_history
user3666197
  • 1
  • 6
  • 50
  • 92
Barnaby Searle
  • 21
  • I don't think Numpy will help you: you need the previous position to calculate the new one. Numpy operates on many variables (e.g., positions) at once, but since you can only calculate them one by one, you can't calculate them all at tones. Perhaps you can try Cython, but be warned that may be thrown into the deeps using that. –  Oct 10 '14 at 11:39
  • If you want to calculate it with your code, go with Everts suggestion. If you are interested in a fast solution of your ode, have a look at [odeint](http://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.odeint.html#scipy.integrate.odeint) and reformulate your problem accordingly. – Bort Oct 10 '14 at 11:52
  • Currently your planet just accelerates in a straight line, so no way the code would detect some kind of "cycle" and stop running `:)` –  Oct 10 '14 at 12:53

2 Answers2

2

It seems like you're stuck in an infinite loop.

abs(x) gives the absolute value of a number which is always positive. This means it will never fall below -0.1 so your while loop will never end.

You either need to change abs(x) or -0.1 to something else.

Zulu
  • 8,765
  • 9
  • 49
  • 56
Anon
  • 21
  • 1
0

One suggestion is to define the "history" arrays before your loop rather than appending to them every iteration. This way, a block of memory is reserved beforehand. This should improve performance. However you would need to work out the size of v_history and x_history in advance.

For example, using 1D arrays:

v_history = np.zeros((k,))
x_history = np.zeros((k,))

where (k,) is the shape of the array.

You would then need to use an index value to store the calculated values in the array

# outside loop
x=0

# inside loop
v_history[x] = v
x +=1

Also you might want to start thinking about broadcasting

Community
  • 1
  • 1
Lee
  • 29,398
  • 28
  • 117
  • 170