Is there a way to define a float array in Python?

Question

For my astronomy homework, I need to simulate the elliptical orbit of a planet around a sun. To do this, I need to use a for loop to repeatedly calculate the motion of the planet. However, every time I try to run the program, I get the following error:

RuntimeWarning: invalid value encountered in power
r=(x**2+y**2)**1.5
Traceback (most recent call last):
File "planetenstelsel3-4.py", line 25, in <module>
ax[i] = a(x[i],y[i])*x[i]
ValueError: cannot convert float NaN to integer

I've done some testing, and I think the problem lies in the fact that the values that are calculated are greater than what fits in an integer, and the arrays are defined as int arrays. So if there was a way do define them as float arrays, maybe it would work. Here is my code:

import numpy as np
import matplotlib.pyplot as plt

dt = 3600 #s
N = 5000
x = np.tile(0, N)
y = np.tile(0, N)
x[0] = 1.496e11 #m
y[0] = 0.0
vx = np.tile(0, N)
vy = np.tile(0, N)
vx[0] = 0.0
vy[0] = 28000 #m/s
ax = np.tile(0, N)
ay = np.tile(0, N)
m1 = 1.988e30 #kg
G = 6.67e-11 #Nm^2kg^-2

def a(x,y):
    r=(x**2+y**2)**1.5
    return (-G*m1)/r

for i in range (0,N):
    r = x[i],y[i]
    ax[i] = a(x[i],y[i])*x[i]
    ay[i] = a(x[i],y[i])*y[i]
    vx[i+1] = vx[i] + ax[i]*dt
    vy[i+1] = vy[i] + ay[i]*dt
    x[i+1] = x[i] + vx[i]*dt 
    y[i+1] = y[i] + vy[i]*dt    
plt.plot(x,y)
plt.show()

The first few lines are just some starting parameters.

Thanks for the help in advance!

Answer 1

When you are doing physics simulations you should definitely use floats for everything. 0 is an integer constant in Python, and thus np.tile creates integer arrays; use 0.0 as the argument to np.tile to do floating point arrays; or preferably use the np.zeros(N) instead:

You can check the datatype of any array variable from its dtype member:

>>> np.tile(0, 10).dtype
dtype('int64')
>>> np.tile(0.0, 10).dtype
dtype('float64')
>>> np.zeros(10)
array([ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.])
>>> np.zeros(10).dtype
dtype('float64')

To get a zeroed array of float32 you'd need to give a float32 as the argument:

>>> np.tile(np.float32(0), 10)
array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0.], dtype=float32)

or, preferably, use zeros with a defined dtype:

>>> np.zeros(10, dtype='float32')
array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0.], dtype=float32)

Answer 2

You need x = np.zeros(N), etc.: this declares the arrays as float arrays.

This is the standard way of putting zeros in an array (np.tile() is convenient for creating a tiling with a fixed array).