Geodesic Trajectories on a Torus

Computational Analysis and Visualization in Differential Geometry

I’ll provide an example problem in differential geometry, solve it using $Python$, and visualize the results with graphs.

I’ll include both the code and mathematical expressions.

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm

# Define the torus parametrization
def torus(u, v, R=2, r=1):
    """
    Parametrize a torus with major radius R and minor radius r
    u, v are the parameters in [0, 2π)
    """
    x = (R + r * np.cos(v)) * np.cos(u)
    y = (R + r * np.cos(v)) * np.sin(u)
    z = r * np.sin(v)
    return x, y, z

# Compute the metric coefficients for the torus
def torus_metric(u, v, R=2, r=1):
    """
    Compute the metric coefficients E, F, G for the torus at point (u, v)
    """
    E = (R + r * np.cos(v))**2  # g_11
    F = 0                       # g_12 = g_21
    G = r**2                    # g_22
    return E, F, G

# Compute the Christoffel symbols for the torus
def christoffel_symbols(u, v, R=2, r=1):
    """
    Compute the Christoffel symbols for the torus at point (u, v)
    Returns a 2x2x2 array where Gamma[i,j,k] = Gamma^i_{jk}
    """
    E, F, G = torus_metric(u, v, R, r)
    
    # Derivatives of metric coefficients
    dE_du = 0
    dE_dv = -2 * (R + r * np.cos(v)) * r * np.sin(v)
    dF_du = dF_dv = 0
    dG_du = 0
    dG_dv = 0
    
    # Initialize Christoffel symbols
    Gamma = np.zeros((2, 2, 2))
    
    # g^{11} = 1/E, g^{12} = g^{21} = 0, g^{22} = 1/G
    g_inv_11 = 1/E
    g_inv_22 = 1/G
    
    # Gamma^1_{11}
    Gamma[0, 0, 0] = 0.5 * g_inv_11 * dE_du
    
    # Gamma^1_{12} = Gamma^1_{21}
    Gamma[0, 0, 1] = Gamma[0, 1, 0] = 0.5 * g_inv_11 * dE_dv
    
    # Gamma^1_{22}
    Gamma[0, 1, 1] = -0.5 * g_inv_11 * dG_du
    
    # Gamma^2_{11}
    Gamma[1, 0, 0] = -0.5 * g_inv_22 * dE_dv
    
    # Gamma^2_{12} = Gamma^2_{21}
    Gamma[1, 0, 1] = Gamma[1, 1, 0] = 0.5 * g_inv_22 * dF_du
    
    # Gamma^2_{22}
    Gamma[1, 1, 1] = 0.5 * g_inv_22 * dG_dv
    
    return Gamma

# Geodesic equation as a first-order system
def geodesic_equation(t, y, R=2, r=1):
    """
    The geodesic equation as a first-order system
    y = [u, v, u_dot, v_dot]
    Returns [u_dot, v_dot, u_ddot, v_ddot]
    """
    u, v, u_dot, v_dot = y
    
    # Compute Christoffel symbols
    Gamma = christoffel_symbols(u, v, R, r)
    
    # Geodesic equations (second derivatives)
    u_ddot = -Gamma[0, 0, 0] * u_dot**2 - 2 * Gamma[0, 0, 1] * u_dot * v_dot - Gamma[0, 1, 1] * v_dot**2
    v_ddot = -Gamma[1, 0, 0] * u_dot**2 - 2 * Gamma[1, 0, 1] * u_dot * v_dot - Gamma[1, 1, 1] * v_dot**2
    
    return [u_dot, v_dot, u_ddot, v_ddot]

# Solve the geodesic equation with given initial conditions
def solve_geodesic(u0, v0, u_dot0, v_dot0, t_span, R=2, r=1):
    """
    Solve the geodesic equation with initial conditions
    u0, v0: initial position
    u_dot0, v_dot0: initial velocity
    t_span: time span [t_start, t_end]
    """
    y0 = [u0, v0, u_dot0, v_dot0]
    
    # Solve the ODE system
    sol = solve_ivp(
        lambda t, y: geodesic_equation(t, y, R, r),
        t_span,
        y0,
        method='RK45',
        rtol=1e-8,
        atol=1e-8,
        dense_output=True
    )
    
    return sol

# Plot the torus and geodesic
def plot_torus_and_geodesic(sol, R=2, r=1, num_points=1000):
    """
    Plot the torus surface and the geodesic curve
    """
    # Create figure
    fig = plt.figure(figsize=(12, 10))
    ax = fig.add_subplot(111, projection='3d')
    
    # Create a mesh grid for the torus
    u = np.linspace(0, 2*np.pi, 100)
    v = np.linspace(0, 2*np.pi, 100)
    u_grid, v_grid = np.meshgrid(u, v)
    
    # Compute the torus coordinates
    x, y, z = torus(u_grid, v_grid, R, r)
    
    # Plot the torus surface with translucent appearance
    ax.plot_surface(x, y, z, cmap=cm.viridis, alpha=0.6, edgecolor='none')
    
    # Extract the solution at evenly spaced points
    t = np.linspace(sol.t[0], sol.t[-1], num_points)
    y_interp = sol.sol(t)
    u_sol, v_sol = y_interp[0], y_interp[1]
    
    # Convert the geodesic curve to Cartesian coordinates
    x_geo, y_geo, z_geo = torus(u_sol, v_sol, R, r)
    
    # Plot the geodesic curve
    ax.plot(x_geo, y_geo, z_geo, 'r-', linewidth=2, label='Geodesic')
    
    # Mark the start and end points
    ax.plot([x_geo[0]], [y_geo[0]], [z_geo[0]], 'go', markersize=8, label='Start')
    ax.plot([x_geo[-1]], [y_geo[-1]], [z_geo[-1]], 'bo', markersize=8, label='End')
    
    # Set labels and title
    ax.set_xlabel('X')
    ax.set_ylabel('Y')
    ax.set_zlabel('Z')
    ax.set_title('Geodesic on a Torus (R={}, r={})'.format(R, r))
    ax.legend()
    
    plt.tight_layout()
    plt.show()
    
    return fig, ax

# Analyze the length and curvature of the geodesic
def analyze_geodesic(sol, R=2, r=1, num_points=1000):
    """
    Analyze properties of the geodesic such as length and curvature
    """
    # Extract the solution at evenly spaced points
    t = np.linspace(sol.t[0], sol.t[-1], num_points)
    y_interp = sol.sol(t)
    u_sol, v_sol = y_interp[0], y_interp[1]
    u_dot_sol, v_dot_sol = y_interp[2], y_interp[3]
    
    # Compute the metric coefficients along the geodesic
    E, F, G = [], [], []
    for i in range(len(t)):
        E_i, F_i, G_i = torus_metric(u_sol[i], v_sol[i], R, r)
        E.append(E_i)
        F.append(F_i)
        G.append(G_i)
    
    E = np.array(E)
    F = np.array(F)
    G = np.array(G)
    
    # Compute the speed along the geodesic
    speed = np.sqrt(E * u_dot_sol**2 + 2*F * u_dot_sol * v_dot_sol + G * v_dot_sol**2)
    
    # Compute the length of the geodesic
    dt = t[1] - t[0]
    length = np.sum(speed * dt)
    
    # Plot the speed along the geodesic
    plt.figure(figsize=(10, 6))
    plt.plot(t, speed)
    plt.title('Speed along the Geodesic')
    plt.xlabel('t')
    plt.ylabel('Speed')
    plt.grid(True)
    plt.show()
    
    print(f"Geodesic Length: {length:.6f}")
    
    return length

# Main execution code
if __name__ == "__main__":
    # Parameters for the torus
    R = 2.0  # Major radius
    r = 1.0  # Minor radius
    
    # Initial conditions for the geodesic
    u0 = 0.0
    v0 = 0.0
    u_dot0 = 1.0
    v_dot0 = 0.5
    
    # Time span for integration
    t_span = [0, 10]
    
    # Solve the geodesic equation
    sol = solve_geodesic(u0, v0, u_dot0, v_dot0, t_span, R, r)
    
    # Plot the results
    plot_torus_and_geodesic(sol, R, r)
    
    # Analyze the geodesic
    geodesic_length = analyze_geodesic(sol, R, r)
    
    # Plot a few different geodesics
    plt.figure(figsize=(12, 10))
    
    initial_conditions = [
        (0.0, 0.0, 1.0, 0.0, 'r', 'Horizontal'),
        (0.0, 0.0, 0.0, 1.0, 'g', 'Vertical'),
        (0.0, 0.0, 1.0, 0.5, 'b', 'Diagonal'),
        (0.0, 0.0, 1.0, 2.0, 'purple', 'Steep Diagonal'),
    ]
    
    ax = plt.axes(projection='3d')
    
    # Plot the torus wireframe
    u = np.linspace(0, 2*np.pi, 30)
    v = np.linspace(0, 2*np.pi, 30)
    u_grid, v_grid = np.meshgrid(u, v)
    x, y, z = torus(u_grid, v_grid, R, r)
    ax.plot_surface(x, y, z, color='lightgray', alpha=0.2)
    
    # Compute and plot different geodesics
    for u0, v0, u_dot0, v_dot0, color, label in initial_conditions:
        sol = solve_geodesic(u0, v0, u_dot0, v_dot0, t_span, R, r)
        
        # Extract the solution
        t = np.linspace(sol.t[0], sol.t[-1], 1000)
        y_interp = sol.sol(t)
        u_sol, v_sol = y_interp[0], y_interp[1]
        
        # Convert to Cartesian coordinates
        x_geo, y_geo, z_geo = torus(u_sol, v_sol, R, r)
        
        # Plot the geodesic
        ax.plot(x_geo, y_geo, z_geo, color=color, linewidth=2, label=label)
    
    ax.set_xlabel('X')
    ax.set_ylabel('Y')
    ax.set_zlabel('Z')
    ax.set_title('Different Geodesics on a Torus (R={}, r={})'.format(R, r))
    ax.legend()
    
    plt.tight_layout()
    plt.show()

Differential Geometry Example: Geodesics on a Torus

In this example, I’ll solve the geodesic equations on a torus surface and visualize the results.

Geodesics are curves that locally minimize distance between points on a manifold, representing the “straightest possible paths” on the surface.

Mathematical Background

Torus Parameterization

A torus with major radius $R$ and minor radius $r$ can be parameterized as:

$$x(u,v) = (R + r\cos v)\cos u$$
$$y(u,v) = (R + r\cos v)\sin u$$
$$z(u,v) = r\sin v$$

where $u, v \in [0, 2\pi)$.

Metric Tensor

The components of the first fundamental form (metric tensor) are:

$$E = g_{11} = (R + r\cos v)^2$$
$$F = g_{12} = g_{21} = 0$$
$$G = g_{22} = r^2$$

Geodesic Equations

The geodesic equations are given by:

$$\frac{d^2u}{dt^2} + \Gamma^1_{11}\left(\frac{du}{dt}\right)^2 + 2\Gamma^1_{12}\frac{du}{dt}\frac{dv}{dt} + \Gamma^1_{22}\left(\frac{dv}{dt}\right)^2 = 0$$

$$\frac{d^2v}{dt^2} + \Gamma^2_{11}\left(\frac{du}{dt}\right)^2 + 2\Gamma^2_{12}\frac{du}{dt}\frac{dv}{dt} + \Gamma^2_{22}\left(\frac{dv}{dt}\right)^2 = 0$$

where $\Gamma^i_{jk}$ are the Christoffel symbols.

Code Explanation

1. Surface Parameterization: I defined the torus using the standard parameterization with major radius $R$ and minor radius $r$.

2. Metric Calculation: I computed the components of the metric tensor ($E$, $F$, $G$).

3. Christoffel Symbols: I calculated the Christoffel symbols from the metric and its derivatives.

4. Geodesic Equations: I formulated the geodesic equations as a first-order ODE system.

5. Numerical Integration: I used solve_ivp from $SciPy$ to numerically integrate the geodesic equations. This converts our second-order differential equations into a system of first-order ODEs.

6. Visualization: I created 3D plots showing the torus surface and the geodesic curves on it.

7. Analysis: I calculated and plotted the speed along the geodesic and computed its total length.

Key Functions in the Code

• torus(u, v, R, r): Computes the Cartesian coordinates (x, y, z) for given parameters (u, v)
• torus_metric(u, v, R, r): Calculates the metric tensor components E, F, G
• christoffel_symbols(u, v, R, r): Computes the Christoffel symbols for the torus
• geodesic_equation(t, y, R, r): Implements the geodesic equations as a first-order system
• solve_geodesic(u0, v0, u_dot0, v_dot0, t_span, R, r): Solves the geodesic initial value problem
• plot_torus_and_geodesic(sol, R, r, num_points): Visualizes the torus and geodesic
• analyze_geodesic(sol, R, r, num_points): Computes geodesic properties like length and speed

Results Interpretation

When you run the code, you’ll see several visualizations:

1. Main Geodesic Visualization:
   
   Shows a single geodesic curve (in red) on the torus surface, with green and blue markers for the start and end points.

2. Speed Profile:
   
   A 2D plot showing how the speed of a particle following the geodesic changes over time.
   For a geodesic with proper parameterization, this speed should be constant, but numerical integration can introduce small variations.

3. Multiple Geodesics:
   
   A comparison of different geodesics with various initial velocities:

   • Horizontal geodesic ($u_dot0 = 1, v_dot0 = 0$)
   • Vertical geodesic ($u_dot0 = 0, v_dot0 = 1$)
   • Diagonal geodesic ($u_dot0 = 1, v_dot0 = 0.5$)
   • Steep diagonal geodesic ($u_dot0 = 1, v_dot0 = 2.0$)

Mathematical Insights

1. Horizontal and Vertical Geodesics: Special cases where the geodesics follow circles on the torus.
   When $u_dot0 = 1$ and $v_dot0 = 0$, the geodesic follows a circle of constant v.
   When $u_dot0 = 0$ and $v_dot0 = 1$, it follows a circle of constant u.

2. Diagonal Geodesics: More complex paths that wind around the torus, potentially never closing back on themselves (creating dense coverage if extended long enough).

3. Geodesic Behavior: Notice how geodesics tend to avoid the inner part of the torus and prefer to travel along the outer part where the curvature is less extreme.
   This is a direct consequence of the geodesic equations trying to minimize the path length.

Applications

This type of analysis has applications in:

• Computer graphics (for path planning on surfaces)
• General relativity (geodesics represent the paths of free-falling particles)
• Robotics (finding optimal paths on curved surfaces)
• Differential geometry (studying the intrinsic properties of surfaces)

You can modify the code to explore different surfaces by changing the parameterization, metric, and Christoffel symbols accordingly.