RuntimeWarning: деление на ноль, встречающееся в true_divide - PullRequest
Новая
       34

RuntimeWarning: деление на ноль, встречающееся в true_divide

1 голос
/ 03 июля 2019

Так что в основном у меня есть значения (координаты), которые меняются со временем.Я определил их в func () и решил их с помощью значений odeint и start.Теперь я пытаюсь проверить, соответствуют ли они условию H = 0 (H определено в коде) в любой момент, в противном случае мои траектории неверны.Однако я не понимаю, откуда возникла ошибка, я не понимаю, почему это не работает.У кого-нибудь есть идея?Если вам нужна дополнительная информация, я программирую легкие орбиты вокруг черной дыры Шварцшильда.(если H ≠ 0, это не вращение фотонов).

Я пытался изменить значения и переместить скобки, но, честно говоря, я действительно не знаю, что еще попробовать 

import numpy as np
from scipy.integrate import odeint
from sympy import *
import sympy as sp



# coordinates who won't change
ptheta = 0
rs = 2
t0 = 0
theta = sp.pi/2

def func(f, T, rs):
       # f
       t, r, theta, phi, pt, pr, ptheta, pphi = f
       # derivatives
       dtdT = (1+(rs/r))*pt # 1/(1-rs/r)
       drdT = (1-(rs/r))*pr
       dthetadT = 0 #ptheta/(r**2)
       dphidT = pphi/((r**2))#*(sp.sin(theta)**2))
       dptdT = 0
       dprdT = ((pphi**2)/(r**3))-((rs/2)/((r-rs)**2))-(((rs/2)*    (pr**2))/(r**2))
       dpthetadT = 0 #(cos(theta)*(ptheta)**2)/((r**2)*sp.sin(phi)**3)
       dpphidT = 0 #pphi/((r**2)*(sp.sin(theta)**2)) #0
       return(dtdT, drdT, dthetadT, dphidT, dptdT, dprdT, dpthetadT, dpphidT)

# initial conditions  
r0 = 12*rs
phi0 = sp.pi
pr0 = 0
pphi0 = 5 
pt0 = sqrt((((1-(rs/r0))**2)*(pr0**2)) + (pphi0**2)*(1-rs/r0)/(r0**2))
f0=[t0, r0, theta, phi0, pt0, pr0, ptheta, pphi0]
T = np.linspace(0, 2250, 9000)
xx=odeint(func, f0, T, args=(rs,))

# calling coordinates from odeint
r00 = xx[:, 1]
theta00 = xx[:, 2]
phi00 = xx[:, 3]
pt00 = xx[:, 4]
pr00 = xx[:, 5]
pphi00 = xx[:, 7]

# hamiltonian
H0 = -((1-rs/r00)**-1)*(pt00**2)/2 + (1-rs/r00)*(pr00**2)/2 +    (pphi00**2)/(2*(r00**2))

print(H0)
if np.amax(H0) < 6e-12:
print("Your results are correct")
else:
print("Your results are wrong")

Я должен получить действительно маленькие значения Гамильтона (x e-17 или меньше), но я продолжаю получать n e-02, который слишком велик, чтобы быть правильным

...