为什么我的光子没有以正确的半径运行？

Question

在 Blender 中，我正在路径追踪一些 2D 光子弯曲，这是由于史瓦西度量定义的非旋转标准黑洞造成的。我根据仿射参数 lambda 的 r（半径）、phi（角度）和 t（时间）设置初始速度，然后根据相应点处的 Christoffel 符号迭代更新时空向量。只需将代码粘贴到 Blender python 脚本中并运行即可。

当从远处投射时，我的光子不会以 3GM/(c^2) 轨道运行。作为参考，内球体是半径为 2GM/(c^2) 的黑洞。在环的右侧，光子进入轨道半径内，并且不会螺旋到事件视界！

然而，当最初从 x = 3GM/(c^2) 投影时，它确实进入了正确半径的轨道！ （最终逃脱，因为它很容易出现轻微的错误）

一切都按比例缩小了 1e6 以在 Blender 中绘制

import bpy
from math import sin, cos, pi, sqrt, atan2
from mathutils import Vector

# 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 = 4000  # Resolution of the curve
#test_point = Vector((0, 4, 0))  # Initial position 
test_point = Vector((3*G*M/(1e6 * c**2), 0, 0))
alpha = 4 * pi/8 # initial projection angle (at x = 0, y = 4 try 2.55 * pi/8)
d_lambda = 1e-4  # Step size for integration

# Initialize variables
r = test_point.length * 1e6  # Convert to meters
p = atan2(test_point.y, test_point.x)  # p = phi = Initial planar angle

# initial radial and angular velocities
v_r = -c * cos(alpha)  # Radial component of velocity
v_p = c * sin(alpha) / r  # Angular component of velocity

# Schwarzschild metric components (2D)
g_tt = -(1 - 2 * G * M / (r * c**2)) # time
g_rr = 1 / (1 - 2 * G * M / (r * c**2)) # radial
g_pp = r**2 # phi

# Initialize velocity in the time direction 
v_t = sqrt(-(g_rr * v_r**2 + g_pp * v_p**2) / g_tt)

# Debug: Initial null geodesic condition
print(g_tt * v_t**2 + g_rr * v_r**2 + g_pp * v_p**2)

# Plot black hole sphere mesh
horizon = 2 * G * M / c**2  # Event horizon radius
bpy.ops.mesh.primitive_uv_sphere_add()
black_hole = bpy.context.object
black_hole.scale = (horizon / 1e6, horizon / 1e6, horizon / 1e6)  # Scale for Blender units

# Create a 2D trajectory curve
bpy.ops.mesh.primitive_circle_add(radius=1, vertices=curve_num + 1)
curve = bpy.context.object
verts = curve.data.vertices

# Delete unnecessary vertex to make the trajectory a 2D curve
bpy.ops.object.mode_set(mode='EDIT')
bpy.ops.mesh.select_all(action='DESELECT')
bpy.ops.object.mode_set(mode='OBJECT')
verts[-1].select = True
bpy.ops.object.mode_set(mode='EDIT')
bpy.ops.mesh.delete(type='VERT')
bpy.ops.object.mode_set(mode='OBJECT')

# Main simulation loop
for i, v in enumerate(verts):
    if i == 0:
        v.co = test_point
    else:
        if r > horizon + 1000:
            # Relevant 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

            # Get accelerations
            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

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

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

            # Recalculate 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

            # Debug: Print null geodesic condition
            print(g_tt * v_t**2 + g_rr * v_r**2 + g_pp * v_p**2)

        # Convert spherical to Cartesian coordinates (2D in equatorial plane)
        x = (r / 1e6) * cos(p)
        y = (r / 1e6) * sin(p)

        # Update vertex position
        v.co = Vector((x, y, 0))

# Finalize the curve
bpy.context.view_layer.objects.active = curve
bpy.ops.object.convert(target='CURVE')
curve = bpy.context.object
curve.data.bevel_depth = 0.01
curve.data.fill_mode = 'FULL'
curve.data.bevel_resolution = 12

我希望 for 循环内的空测地线打印语句为 0 或相对较小的整数。例如，第一个空测地线打印语句可能是 0、16、-8 等，因为大浮点加法（e17 幅度）存在少量不精确性。

我尝试通过将“else”语句替换为“elif j == 1”来进行调试，并查看下一次迭代，可以看到空测地线打印出更大的浮点数。

我相信找出零测地线误差将揭示为什么我的轨迹不正确。我已经多次检查了我的 Christoffel 符号和实现，但没有发现任何明显的东西。

Answer 1

您在加速度分量的表达式中缺少因子 2

dv_p

。将该行更改为

dv_p = - 2 * Γ_p_rp * v_r * v_p

这是因为

γ_p_pr * v_p * v_r = γ_p_rp * v_r * v_p

并且您还需要在求和中包含该术语。