处理轨迹奇异性：删除光子

Question

我正在编码一个黑洞（实际上是光子围绕黑洞运行），我需要处理半径值小于极限距离的异常

我试过使用if和while True

def Hamiltonian(r, pt, pr, pphi):
    H = (-((1-rs/r)**-1)*(pt**2)/2 + (1-rs/r)*(pr**2)/2 + (pphi**2)/(2*         (r**2)))
    if np.amax(H) < 10e-08:
        print("Your results are correct")
    else:
        print("Your results are wrong")
    return(H)    

def singularity(H, r):
    if (r).any < 1.5*rs:
        print(H)
    else:
        print("0")
        return(H, r)

print(Hamiltonian(r00, pt00, pr00, pphi00))   

我想处理r<；1.5*rs的情况，这样我就不再有除法错误信息和奇怪的轨道了。到目前为止，我的代码没有改变任何问题，我仍然得到以下错误消息：

"RuntimeWarning: divide by zero encountered in true_divide
  H = (-((1-rs/r)**-1)*(pt**2)/2 + (1-rs/r)*(pr**2)/2 + (pphi**2)/(2*(r**2)))"

我的轨道是完全错误的（例如，我的光子应该在黑洞中右转，但由于在r<；1.5*rs处有一个奇点，它们进入然后又离开，以此类推）

我想删除导致问题的光子，但我不知道如何删除，有人能帮我吗？你知道吗

Answer 1

我想你提到我们不知道奇点到底发生了什么。这里提供的任何答案很可能都不准确，但我们假设您知道接近0的邻居的行为/动态。我不知道你怎么称呼哈密顿函数，但我可以想象你是用两种方法中的一种。你知道吗

1. 您有一个预定义的循环，该循环遍历传递到函数中的每个值，并为该预定义循环中的每个值输出结果H。

2. 您正在传递相同长度的向量，并在H函数中逐个元素地进行数学运算。

在第一种情况下，您可以进行预检查并为接近0的行为编写一个新函数，如果您在接近0的邻域中，则可以调用该新函数。你也可以检查哈密顿函数本身，如果你在接近0的邻域中，调用新函数。后一种方法是我更喜欢的方法，因为它保持了前端（我对它最好的描述）相当干净，并将逻辑/数学封装在哈密顿函数中。在你的评论中，你说你想删除某个半径范围内的光子，这将非常容易用这种方法来完成，只要终止for循环，在那个时间步中断并绘制你在那个时间步之前拥有的东西。你知道吗

在第二种情况下，如果向量属于奇点邻域，则必须通过检查整个向量来手动构建新向量。这将是一个更困难的一点，将取决于您的输入形状和做元素的元素对元素的数学是更难调试在我看来。你知道吗