Here is a thought; Instead of using F=ma, use the equations of motion from GR:

================================================================================ SCHWARZSCHILD METRIC AND GEODESIC EQUATIONS OF MOTION (SUMMARY) ================================================================================

I. THE SCHWARZSCHILD METRIC (g_uv)

The spacetime geometry is defined by the *line element*, ds^2, which relates coordinate changes (dt, dr, d(phi), etc.) to physical distance or proper time: ds^2 = g_uv * dx^u * dx^v

For the Schwarzschild vacuum solution, the line element in the equatorial plane (theta = pi/2) is: ds^2 = -(1 - r_s / r) * c^2 * dt^2 + (1 - r_s / r)^(-1) * dr^2 + r^2 * d(phi)^2

The corresponding non-zero metric components (g_uv) are: g_tt = -(1 - r_s / r) * c^2 g_rr = 1 / (1 - r_s / r) g_phiphi = r^2

Where: r_s = 2 G * M / c^2 (Schwarzschild Radius)

The Lagrangian L for the geodesic path is constructed directly from the metric: L = (1/2) * [ g_tt * (dt/d(lambda))^2 + g_rr * (dr/d(lambda))^2 + g_phiphi (d(phi)/d(lambda))^2 ]

--------------------------------------------------------------------------------

II. CONSERVATION LAWS (FROM EULER-LAGRANGE EQUATIONS)

A. TIME EOM (Conserved Energy E) Since the metric is time-independent, the quantity conjugate to t is conserved: *Specific Energy (E)*.

EQUATION (1): Time Evolution d(t)/d(lambda) = E / ( c^2 * (1 - r_s / r) )

B. PHI EOM (Conserved Angular Momentum L_z) Since the metric is symmetric with respect to phi, the quantity conjugate to phi is conserved: *Specific Angular Momentum (L_z)*.

EQUATION (2): Angular Evolution d(phi)/d(lambda) = L_z / r^2

--------------------------------------------------------------------------------

III. RADIAL EQUATION OF MOTION (FROM THE METRIC CONSTRAINT)

The radial EOM is derived by imposing the metric normalization condition (g_uv * u^u * u^v = epsilon).

A. MASSIVE PARTICLES (Mass m > 0) The proper time (tau) is the affine parameter (lambda=tau), and the normalization is epsilon = c^2. The final EOM is: (dr/d(tau))^2 = E^2/c^2 - V_eff^2

EQUATION (3M): Radial EOM (Massive) (dr/d(tau))^2 = E^2/c^2 - c^2 * (1 - r_s/r) * ( 1 + L_z^2 / (c^2 * r^2) )

B. MASSLESS PARTICLES (Mass m = 0) The normalization is epsilon = 0. The final EOM is: (dr/d(lambda))^2 = E^2 - V_eff^2

EQUATION (3P): Radial EOM (Massless / Photon) (dr/d(lambda))^2 = E^2 - (1 - r_s/r) * L_z^2 / r^2

--------------------------------------------------------------------------------

IV. SUMMARY OF GEODESIC EQUATIONS OF MOTION (EOM)

The motion of any particle (massive or massless) in the Schwarzschild spacetime is determined by the following three coupled first-order differential equations:

A. TIME EVOLUTION: d(t)/d(lambda) = E / ( c^2 * (1 - r_s / r) )

B. ANGULAR EVOLUTION: d(phi)/d(lambda) = L_z / r^2

C. RADIAL EVOLUTION (Specific): 1. Massive Particle (using d(tau)): (dr/d(tau))^2 = E^2/c^2 - c^2 * (1 - r_s/r) * ( 1 + L_z^2 / (c^2 * r^2) )

2. Massless Particle (using d(lambda)): (dr/d(lambda))^2 = E^2 - (1 - r_s/r) * L_z^2 / r^2

================================================================================

This also holds for a non-rotating black hole.

N-Body Simulator – Interactive 3 Body Problem and Gravitational Physics