Почему мой фотон не движется по орбите с правильным радиусом? ⇐ Python

[Anonymous]

В Blender я отслеживаю путь некоторых двумерных фотонов, изгибающихся из-за невращающейся стандартной черной дыры, определенной метрикой Шварцшильда. Я устанавливаю свои начальные скорости в терминах r (радиус), фи (угол) и t (время) относительно аффинного параметра лямбда, а затем итеративно обновляю вектор пространства-времени на основе символов Кристоффеля в соответствующей точке. Просто вставьте код в скрипт Python Blender и запустите его.
При проецировании издалека мой фотон не вращается по орбите 3GM/(c^2). Для справки: внутренняя сфера — это черная дыра радиуса 2GM/(c^2). В правой части петли фотон входит в радиус орбиты и не движется по спирали к горизонту событий!


Однако, если изначально проецировать из x = 3GM/(c^2), это действительно происходит на орбиту правильного радиуса! (Со временем ускользает, так как уязвим к небольшим неточностям)


Все уменьшено на 1e6 для построения в Blender

Код: Выделить всё

import numpy as np
import matplotlib.pyplot as plt
from math import sin, cos, sqrt, atan2, pi

# Constants
M = 9e32  # Mass of the black hole in kg
c = 299792458  # Speed of light in m/s
G = 6.6743e-11  # Gravitational constant in N m^2/kg^2

# Simulation parameters
curve_num = 2000  # Resolution of the curve
d_lambda = 1e-4  # Step size for integration

# Black hole parameters
horizon = 2 * G * M / c**2  # Event horizon radius

# Test point parameters
num = 10 # number of photons
x_i = 5
y_i = 4.7
y_f = 5.5

# Initial velocity in Cartesian coordinates
v_x = -c  # Speed of light in negative x-direction
v_y = 0   # No velocity in the y-direction

# Prepare for matplotlib plotting
trajectories = []

for i in range(num):
trajectory = []

# Define initial test point
if num > 1:
test_point = np.array([x_i, y_i + i / (num - 1) * (y_f - y_i), 0])
else:
test_point = np.array([x_i, y_i, 0])

# Convert to polar coordinates
r = np.linalg.norm(test_point) * 1e6  # Radius in meters
p = atan2(test_point[1], test_point[0])  # Initial planar angle

# Metric coefficients
g_tt = -(1 - 2 * G * M / (r * c**2))
g_rr = 1 / (1 - 2 * G * M / (r * c**2))
g_pp = r**2

# Initial velocities
v_r = (v_x * cos(p) + v_y * sin(p))  # Radial velocity
v_p = (-v_x * sin(p) + v_y * cos(p)) / r  # Angular velocity
v_t = sqrt(-(g_rr * v_r**2 + g_pp * v_p**2) / g_tt)

# Print null geodesic condition for initial conditions
null_condition_initial = g_tt * v_t**2 + g_rr * v_r**2 + g_pp * v_p**2
print(f"Initial null geodesic condition for photon point {i + 1}: {null_condition_initial:.6e}")

# Integrate geodesics
for j in range(curve_num):
if r >  horizon + 1000:
# Christoffel symbols
Γ_r_tt = (G * M / (r**2 * c**2)) * (1 - 2 * G * M / (r * c**2))
Γ_r_rr = -(G * M / (r**2 * c**2)) / (1 - 2 * G * M / (r * c**2))
Γ_r_pp = -r * (1 - 2 * G * M / (r * c**2))
Γ_p_rp = 1 / r
Γ_t_rt = (G * M / (r**2 * c**2))
Γ_t_rt = (G * M / (r**2 * c**2)) / (1 - 2 * G * M / (r * c**2))

# Update change in velocities
dv_r = -Γ_r_tt * v_t**2 - Γ_r_rr * v_r**2 - Γ_r_pp * v_p**2
dv_p = -Γ_p_rp * v_r * v_p
dv_t = -Γ_t_rt * v_r * v_t

# Update velocities
v_r += dv_r * d_lambda
v_p += dv_p * d_lambda
v_t += dv_t * d_lambda

# Update positions
r += v_r * d_lambda
p += v_p * d_lambda

# Check null geodesic condition at each step
null_condition = g_tt * v_t**2 + g_rr * v_r**2 + g_pp * v_p**2
print(f"Step {j + 1}, Photon point {i + 1}, Null condition: {null_condition:.6e}")

# Store Cartesian coordinates
x = (r / 1e6) * cos(p)
y = (r / 1e6) * sin(p)
trajectory.append((x, y))
else:
break

trajectories.append(trajectory)

# Plot using matplotlib
plt.figure(figsize=(8, 8))

# Plot each trajectory
for i, trajectory in enumerate(trajectories):
trajectory = np.array(trajectory)
plt.plot(trajectory[:, 0], trajectory[:, 1], label=f"Photon Trajectory {i + 1}")

# Plot the event horizon
circle = plt.Circle((0, 0), horizon / 1e6, color='black', fill=True, label="Event Horizon")
plt.gca().add_artist(circle)

# Determine plot limits based on trajectory data
all_x = [point[0] for traj in trajectories for point in traj]
all_y = [point[1] for traj in trajectories for point in traj]
max_extent = max(max(map(abs, all_x)), max(map(abs, all_y)))
plt.xlim(-max_extent, max_extent)
plt.ylim(-max_extent, max_extent)

# Configure plot
plt.title("Photon Trajectories Around a Black Hole")
plt.xlabel("x (scaled by 1e6 meters)")
plt.ylabel("y (scaled by 1e6 meters)")
plt.axhline(0, color="gray", linewidth=0.5)
plt.axvline(0, color="gray", linewidth=0.5)
plt.legend()
plt.axis('equal')
plt.grid()
plt.show()

Я ожидаю, что мой оператор нулевой геодезической печати внутри цикла for будет либо 0, либо относительно небольшим целым числом. Например, первый оператор нулевой геодезической печати может иметь значение 0, 16, -8 и т. д. из-за небольшой неточности при большом добавлении чисел с плавающей запятой (величины e17).
У меня есть попытался выполнить отладку, заменив оператор «else» на «elif j == 1» и, посмотрев на следующую итерацию, можно увидеть, что нулевая геодезическая печатает гораздо большее число с плавающей запятой.
Я верю, что подумаю Нулевая геодезическая ошибка покажет, почему мои траектории неверны. Я несколько раз проверял символы и реализацию Кристоффеля и не нашел ничего очевидного.

С пропусками:




Без пропуски:



Подробнее здесь: https://stackoverflow.com/questions/793 ... ect-radius