Как избавиться от ZeroDivisionError: деление на ноль

Question

Я продолжаю сталкиваться с ошибкой деления числа на ноль и не могу понять, почему я это получаю Однако, когда я запускаю код, изначально предоставленный мне (написанный и запущенный в matlab), ошибки не возникают.

Код

import numpy as np
import matplotlib.pyplot as plt
from astropy import constants as const


#Part 1: Exploring Rotation Curves


M = 10**42 #Approximate mass of the Milky Way (kg)
G = const.G #Universal gravitational constant (m^3 kg^-1 s^-2)
r = np.linspace(0, 3e20) #Radii (m)
rkpc = r*(3.24e-20) #Radii (kpc)
plt.figure(1)
plt.title('Rotation Curves for Three Mass Distributions') 
v1 = np.sqrt(G * M / r) # Orbital velocity in system with central mass (m/s)
M_prop = np.linspace(0, M) # Array of masses increasing proportionally with radius
v2 = np.sqrt(G * M_prop / r) 
M_dens = (M * (r / (max(r)))**3) 
v3 = np.sqrt((G * M_dens) / r)
plt.plot(rkpc, v1/1000, 'b', label = 'Constant M_{r}')
plt.plot(rkpc, v2/1000, 'k', label = 'M_{r} \propto r')
plt.plot(rkpc, v3/1000, 'r', label = 'M_{r} \propto r^{3}')

Я знаю, что ошибка происходит из-за двух следующих строк

M_dens = (M * (r / (max(r)))**3) 
v3 = np.sqrt((G * M_dens) / r)

Я предполагаю, что это происходит из-за максимума (r), но сможет ли кто-то пролить больше света на то, почему это происходит? Потенциально исправление?

Answer 1

Я не уверен, как Matlab обрабатывает деление на ноль, но возможно изменить поведение numpy, используя np.errstate.

a = np.arange(-5, 5.)
b = np.arange(-2, 8.)
with np.errstate(divide='ignore'):
     res0 = a / b
     res1 = b / a 
print(res0, '\n', res1)
# [ 2.5    4.   -inf -2.   -0.5   0.  0.25   0.4  0.5  0.57142857] 
# [ 0.4   0.25 -0.   -0.5  -2.   inf  4.    2.5   2.    1.75]

В качестве альтернативы создайте функцию, которая может установить inf, -inf приводит к полезному значению по умолчанию.

def do_div( a,b, def_val=np.inf):
    with np.errstate(divide='ignore'):
        res = a / b                   
        res[ ~np.isfinite(res) ] = def_val
    return res

print( do_div( a, b, 100 ))
# [  2.5   4.   100.  -2.    -0.5   0.     0.25   0.4  0.5  0.57142857]
print( do_div( b, a, 100 ))
# [  0.4  0.25  -0.  -0.5  -2. 100.  4.  2.5  2.  1.75]

Установка состояния деления для деления на «игнорировать» подавляет предупреждение. Numpy возвращает плюс или минус бесконечность для деления на ноль. Функция do_div устанавливает любые значения бесконечности по умолчанию. В моей работе это чаще всего ноль. Я использовал 100 здесь, так что это легко увидеть. Matlab, вероятно, делает нечто подобное, возвращая бесконечность или альтернативное значение по умолчанию и не выдавая ошибку или предупреждение.

Answer 2

Извините, если это не сработает, я немного груб с математическими командами, подобными этим.

В этой строке:

r = np.linspace(0, 3e20)

r начнется как 0. Позже в эта строка:

v3 = np.sqrt((G * M_dens) / r)

вы делите на r, что равно 0. Все, что разделено на 0, не определено, поэтому Python это не нравится и выдает ошибку.