Dealing with extremely small/large values in odeint

Question

I am trying to solve this set of coupled differential equations using scipy.integrate.odeint.

I wrote this Python code:

def odes(x, r):
    #constants
    alpha = 1.473
    beta = 52.46
    gamma = 1.33 
    
    #Assign each ode to a vector element
    p = x[0]
    m = x[1]
    
    #Define each ode
    dp_dt = -alpha*(x[0]**(1/gamma))*x[1]/(r**2)
    dm_dt = beta*(r**2)*(x[0]**(1/gamma))
    
    #Return list
    return [dp_dt, dm_dt]

#initial conditions(Initial p_bar, initial m_bar)
init = [1e-16, 0]

#declare a r vector
r = np.linspace(0, 15000, 1000000)

#Solve the odes
x = odeint(odes, init, r)

print(x)
plt.semilogy(r, x[:,0])
plt.semilogy(r, x[:,1])
plt.show()

However, Jupyter threw this RuntimeWarning:

C:\conda_temp\ipykernel_14612\3000110852.py:12: RuntimeWarning: invalid value encountered in double_scalars
  dp_dt = -alpha*(x[0]**(1/gamma))*x[1]/(r**2)

This is what the x array looks like:

[[1.e-16 0.e+00]
 [   nan    nan]
 [0.e+00 0.e+00]
 ...
 [0.e+00 0.e+00]
 [0.e+00 0.e+00]
 [0.e+00 0.e+00]]

From what I gathered, it has something to do with how Python/odeint handles extremely small values. However, I have no idea how to get around this problem. I have tried using scipy.special.logsumexp (incorrectly?) to no avail.

Answer 1

Your issue isn't the size of the numbers, your issue is that your first r value is 0, which means you are dividing by zero when you compute dp_dt. Instead, set the first r value to something small but non-zero, e.g.

r = np.linspace(1e-6, 15000, 1000000)