Using matplotlib, I am path tracing some 2D photons bending due to a non-rotating standard black hole defined by Schwarzschild metric. I set my initial velocities in terms of r (radius), phi (angle), and t (time) with respect to the affine parameter lambda and then iteratively update the space-time vector based on the Christoffel symbols at the respective point. The main concerns are is why the null geodesic statements don't print something closer to zero. Everything scaled down by 1e6 to plot 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 photon_sphere = 3 * G * M / c**2 # Photon sphere radius # Test point parameters num = 30 # Number of photons x_i = 7.5 y_i = 2.5 y_f = 4 # 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]) # linear interpolation 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) # Check the null geodesic condition for the initial point #print(g_tt * v_t**2 + g_rr * v_r**2 + g_pp * v_p**2) # Integrate geodesics for j in range(curve_num): if r > horizon + 10000: # Precompute common terms term1 = G * M / (r**2 * c**2) # Common term: GM / r^2c^2 term2 = 1 - 2 * G * M / (r * c**2) # Common term: 1 - 2GM / rc^2 # Christoffel symbols using common terms Γ_r_tt = term1 * term2 Γ_r_rr = -term1 / term2 Γ_r_pp = -r * term2 Γ_p_rp = 1 / r #Γ_t_rt = term1 / term2 # ignoring time marching # Update change in velocities dv_r = -Γ_r_tt * v_t**2 - Γ_r_rr * v_r**2 - Γ_r_pp * v_p**2 dv_p = -2 * Γ_p_rp * v_r * v_p #dv_t = -2 * Γ_t_rt * v_r * v_t # ignoring time marching # Update velocities v_r += dv_r * d_lambda v_p += dv_p * d_lambda #v_t += dv_t * d_lambda # ignoring time marching # Update positions r += v_r * d_lambda p += v_p * d_lambda # Metric tensor components (update for new r) g_tt = -term2 g_rr = 1 / term2 g_pp = r**2 # Recalculate v_t from the metric components v_t = sqrt(-(g_rr * v_r**2 + g_pp * v_p**2) / g_tt) # Check the null geodesic condition at each step #print(g_tt * v_t**2 + g_rr * v_r**2 + g_pp * v_p**2) # Store Cartesian coordinates x = (r / 1e6) * cos(p) y = (r / 1e6) * sin(p) # Only store points within the -10 to 10 range if -10 <= x <= 10 and -10 <= y <= 10: trajectory.append((x, y)) else: break trajectories.append(trajectory) # Plot using matplotlib plt.figure(figsize=(8, 8)) # Plot each trajectory for trajectory in trajectories: trajectory = np.array(trajectory) if len(trajectory) > 0: # Plot only if there are points plt.plot(trajectory[:, 0], trajectory[:, 1]) # Plot the event horizon circle = plt.Circle((0, 0), horizon / 1e6, color='black', fill=True, label="Event Horizon") plt.gca().add_artist(circle) # Plot the photon sphere photon_sphere_circle = plt.Circle((0, 0), photon_sphere / 1e6, color='red', fill=False, linestyle='--', linewidth=1.5, label="Photon Sphere") plt.gca().add_artist(photon_sphere_circle) # Set plot limits explicitly plt.xlim(-10, 10) plt.ylim(-10, 10) # 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.axis('equal') plt.grid() plt.legend() plt.show()  I expect my null geodesic print statement inside the for loop to be either 0 or a relatively small integer. For example, the first null geodesic print statement may be 0, 16, -8, etc. due to a small amount of imprecision with the large float addition (e17 magnitudes). I have tried to debug by replacing my "else" statement with "elif j == 1" and looking at the very next iteration, it can be seen the null geodesic prints a much larger float. I believe figuring out the null geodesic error will reveal if my trajectories are incorrect.  
You are missing a factor of 2 in the expression for the acceleration component dv_p. Change that line to dv_p = - 2 * Γ_p_rp * v_r * v_p  This is because Γ_p_pr = Γ_p_rp and you need to include that term as well in the summation. Neil Butcher's post also picks out another non-zero Christoffel symbol Γ_t_rt that I missed and would be necessary to get v_t right and then affect the trajectories. BUT ... I believe the metric check is a red herring - you initialise v_t at the start with that metric and I think you can (and should) do the same again each time the metric components are recalculated from a new r. That is then bound to keep you on a geodesic. For completeness, I believe that you are solving systems of the form  or, in terms of your quasi-velocities,  If you do have to print the null-geodesic error then I think it should be normalised. You do have to see it in the context of c2, which is very large. Anyway, making the factor-of-2 change does reduce this error by 10 orders of magnitude. You will have to tell me the outcome of the plotting. I don't have, or know how to use, Blender. As you define a variable horizon it would be useful to use that (or the more common rs), rather than repeating 2 * G * M / c**2 In summary:  put the extra factor of 2 in the expression for dv_p update v_t from the metric coefficients every time you recalculate the metric from a new r. 