riemannian

riemannian.py — Unified 1D/2D Riemannian geometry toolkit

Overview

The riemannian module provides a unified framework for working with Riemannian manifolds in one and two dimensions. A single class Metric dispatches all computations automatically based on the dimension, making it easy to switch between 1D curves and 2D surfaces without changing the calling interface.

Key features include:

• Symbolic construction of the metric tensor from an explicit expression (1D scalar, 2D matrix) or by extraction from a Hamiltonian kinetic energy.

• Automatic computation of Christoffel symbols (1D and 2D).

• Geodesic integration with multiple numerical schemes (RK4, adaptive, symplectic/Verlet, Hamiltonian flow via the companion symplectic module).

• Curvature: Riemann tensor, Ricci tensor, Gaussian curvature, scalar curvature.

• Laplace–Beltrami operator: full symbol (principal + subprincipal parts) ready for microlocal analysis.

• Riemannian volume (arc length in 1D) via symbolic or numerical integration.

• 2D only: Exponential map, geodesic distance (shooting or optimisation), Jacobi equation solver (geodesic deviation), Hodge star operator, de Rham Laplacian on 0‑ and 1‑forms, numerical verification of the Gauss–Bonnet theorem.

• Rich visualisation suite: geodesic trajectories, curvature maps (Gaussian/Ricci), metric components.

Mathematical background

A Riemannian metric g on an n-dimensional manifold assigns an inner product to each tangent space. In local coordinates (x¹,…,xⁿ) the metric is written as

ds² = gᵢⱼ(x) dxⁱ dxʲ

and its inverse is denoted gⁱʲ. The Christoffel symbols are derived from the metric:

Γⁱⱼₖ = ½ gⁱˡ (∂ⱼ gₖₗ + ∂ₖ gⱼₗ − ∂ₗ gⱼₖ)

and determine the geodesic equation

ẍⁱ + Γⁱⱼₖ ẋʲ ẋᵏ = 0.

For a 1D metric g₁₁(x) the geodesic equation simplifies to

ẍ + Γ¹₁₁ ẋ² = 0, Γ¹₁₁ = ½ (log g₁₁)′.

Curvature is encoded in the Riemann tensor Rⁱⱼₖₗ, the Ricci tensor Rᵢⱼ = Rᵏᵢₖⱼ, and the scalar curvature R = gⁱʲ Rᵢⱼ. For a 2D surface the Gaussian curvature K satisfies R₁₂₁₂ = K |g| and is the only independent component.

The Laplace–Beltrami operator acting on functions is

Δ = |g|^{-½} ∂ᵢ ( |g|^{½} gⁱʲ ∂ⱼ ),

and its principal symbol is gⁱʲ ξᵢ ξⱼ. The subprincipal symbol encodes the lower‑order terms.

The module also implements the Hodge star on differential forms and the de Rham Laplacian Δ = dδ + δd for 0‑ and 1‑forms in 2D.

References

[1]

do Carmo, M. P. Riemannian Geometry, Birkhäuser, 1992.

[2]

Jost, J. Riemannian Geometry and Geometric Analysis, Springer, 2011 (6th ed.).

[3]

Lee, J. M. Riemannian Manifolds: An Introduction to Curvature, Springer, 1997.

[4]

Petersen, P. Riemannian Geometry, Springer, 2016 (3rd ed.).

[5]

Frankel, T. The Geometry of Physics, Cambridge University Press, 2011 (3rd ed.).

class riemannian.Metric(g_input, coords)

Bases: object

Riemannian metric on a 1D or 2D manifold.

The dimension is inferred from the inputs:

• 1D: g_input is a scalar sympy expression; coords is a 1-tuple (x,) or a single symbol.

• 2D: g_input is a 2×2 sympy Matrix (or nested list); coords is a 2-tuple (x, y).

Parameters:

• g_input (sympy expression | sympy Matrix | list) – Metric tensor. Scalar for 1D, 2×2 matrix for 2D.

• coords (tuple of sympy symbols) – Coordinate variables. Length must match dimension.

dim

Manifold dimension (1 or 2).

Type:

int

coords

Coordinate symbols.

Type:

tuple

christoffel_sym

Symbolic Christoffel symbol(s).

Type:

sympy expression | dict

christoffel_func

Numerical Christoffel symbol(s).

Type:

callable | dict

Examples

>>> x = symbols('x', real=True, positive=True)
>>> m = Metric(x**2, (x,))            # 1D cone-like metric
>>> m.gauss_curvature()               # 0 for 1D
0

>>> x, y = symbols('x y', real=True)
>>> g = Matrix([[1, 0], [0, sin(x)**2]])
>>> m = Metric(g, (x, y))             # unit sphere metric
>>> m.gauss_curvature()

arc_length(x_min, x_max, method='numerical')

Arc length between two points (1D only).

Equivalent to riemannian_volume((x_min, x_max), method).

eval(*point)

Evaluate metric quantities numerically at a point.

Parameters:

*point (float or ndarray) – Coordinate values. Pass one argument for 1D, two for 2D.

Returns:

Keys: ‘g’, ‘g_inv’, ‘sqrt_det’, and ‘christoffel’.

Return type:

dict

classmethod from_hamiltonian(H_expr, coords, momenta)

Extract metric from Hamiltonian kinetic term.

For H = ½ g^{ij} pᵢ pⱼ + V, recover gᵢⱼ via the momentum Hessian.

Parameters:

• H_expr (sympy expression) – Full Hamiltonian H(coords, momenta).

• coords (tuple) – Position symbols, length 1 or 2.

• momenta (tuple) – Momentum symbols, same length as coords.

Returns:

With dimension matching len(coords).

Return type:

Metric

Examples

>>> x, p = symbols('x p', real=True)
>>> H = p**2 / (2*x**2)
>>> m = Metric.from_hamiltonian(H, (x,), (p,))
>>> m.dim
1

gauss_curvature()

Gaussian curvature K.

Returns 0 for 1D (intrinsic curvature vanishes on a curve). For 2D, computes K = R₁₂₁₂ / |g|.

Return type:

sympy expression

laplace_beltrami_symbol()

Symbol of the Laplace-Beltrami operator.

Works for both 1D and 2D.

Returns:

Keys: ‘principal’, ‘subprincipal’, ‘full’.

Return type:

dict

ricci_scalar()

Scalar curvature R.

For 1D returns 0; for 2D computes R = gⁱʲ Rᵢⱼ.

ricci_tensor()

Ricci tensor Rᵢⱼ (2D only).

riemann_tensor()

Riemann curvature tensor Rⁱⱼₖₗ (2D only).

Returns:

Nested dict R[i][j][k][l].

Return type:

dict

Raises:

NotImplementedError – For 1D manifolds (tensor is identically zero).

riemannian_volume(domain, method='numerical')

Riemannian volume of a domain.

Parameters:

• domain (varies by dimension) –

  • 1D: (x_min, x_max)

  • 2D: ((x_min, x_max), (y_min, y_max))

• method ({'symbolic', 'numerical'})

Return type:

float or sympy expression

riemannian.christoffel(metric)

Return the numerical Christoffel symbol(s) of the metric.

Parameters:

metric (Metric)

Returns:

• For dim=1 (callable Γ(x).)

• For dim=2 (nested dict Gamma[i][j][k] of callables.)

riemannian.de_rham_laplacian(metric, form_degree)

Hodge Laplacian Δ = dδ + δd on k-forms (2D only).

Parameters:

• metric (Metric (dim=2))

• form_degree ({0, 1})

Returns:

Symbol dict.

Return type:

dict

riemannian.distance(metric, p, q, method='shooting', max_iter=50, tol=1e-06)

Geodesic distance between two points on a 2D manifold.

Parameters:

• metric (Metric (dim=2))

• p (tuple)

• q (tuple)

• method ({'shooting', 'optimize'})

Return type:

float

riemannian.exponential_map(metric, p, v, t=1.0, method='rk45')

Exponential map exp_p(tv) on a 2D manifold.

Parameters:

• metric (Metric (dim=2))

• p (tuple) – Base point (x₀, y₀).

• v (tuple) – Initial tangent vector (vₓ, vᵧ).

• t (float) – Parameter value.

• method (str)

Returns:

End point (x(t), y(t)).

Return type:

tuple

riemannian.geodesic_hamiltonian_flow(metric, p0, v0, tspan, method='verlet', n_steps=1000)

Integrate geodesic flow in Hamiltonian formulation using symplectic.hamiltonian_flow.

H = ½ gⁱʲ pᵢ pⱼ (1D or 2D).

Parameters:

• metric (Metric)

• p0 (float (1D) | tuple (2D)) – Initial position.

• v0 (float (1D) | tuple (2D)) – Initial velocity (converted to momentum internally).

• tspan (tuple (t_start, t_end))

• method (str) – For 1D/2D: ‘verlet’, ‘stormer’, ‘symplectic_euler’, ‘rk45’. (‘stormer’ is aliased to ‘verlet’, ‘symplectic_euler’ to ‘symplectic’.)

• n_steps (int) – Number of time steps.

Returns:

Phase‑space trajectory with keys: 1D: ‘t’, ‘x’, ‘v’, ‘p’, ‘energy’ 2D: ‘t’, ‘x’, ‘y’, ‘vx’, ‘vy’, ‘px’, ‘py’, ‘energy’

Return type:

dict

riemannian.geodesic_solver(metric, p0, v0, tspan, method='rk4', n_steps=1000, reparametrize=False)

Integrate the geodesic equation.

Dispatches to the appropriate 1D or 2D integrator.

Parameters:

• metric (Metric)

• p0 (float (1D) or tuple (x₀, y₀) (2D)) – Initial position.

• v0 (float (1D) or tuple (vₓ₀, vᵧ₀) (2D)) – Initial velocity / tangent vector.

• tspan (tuple) – (t_start, t_end).

• method (str) – ‘1D: rk4’ | ‘adaptive’ | ‘symplectic’; ‘2D: rk45’ | ‘rk4’ | ‘symplectic’ | ‘verlet’.

• n_steps (int)

• reparametrize (bool) – If True (2D only), also compute arc-length parameter.

Returns:

Trajectory. Keys depend on dimension: 1D → ‘t’, ‘x’, ‘v’ (and ‘p’ for symplectic). 2D → ‘t’, ‘x’, ‘y’, ‘vx’, ‘vy’ (and more for Hamiltonian methods).

Return type:

dict

riemannian.hodge_star(metric, form_degree)

Hodge star operator on differential forms (2D only).

Parameters:

• metric (Metric (dim=2))

• form_degree ({0, 1, 2})

Return type:

callable

riemannian.jacobi_equation_solver(metric, geodesic, initial_variation, tspan, n_steps=1000)

Solve the Jacobi equation (geodesic deviation) on a 2D manifold.

D²J/dt² + R(J, γ̇)γ̇ = 0

Parameters:

• metric (Metric (dim=2))

• geodesic (dict) – Output of geodesic_solver.

• initial_variation (dict) – ‘J0’: (Jx₀, Jy₀), ‘DJ0’: (DJx₀, DJy₀).

• tspan (tuple)

• n_steps (int)

Returns:

‘J_x’, ‘J_y’, ‘DJ_x’, ‘DJ_y’ arrays.

Return type:

dict

riemannian.laplace_beltrami(metric)

Return the Laplace-Beltrami symbol dict for a 1D or 2D metric.

Shortcut for metric.laplace_beltrami_symbol().

riemannian.sturm_liouville_reduce(metric, potential_expr=None)

Reduce the Laplace-Beltrami eigenvalue problem to Sturm-Liouville form.

-Δg u + V u = λ u → -(p u’)’ + q u = λ w u

where p = √g g¹¹, w = √g, q = V √g.

Parameters:

• metric (Metric (dim=1))

• potential_expr (sympy expression, optional) – Potential V(x).

Returns:

Keys: ‘p’, ‘q’, ‘w’ (sympy), ‘p_func’, ‘q_func’, ‘w_func’.

Return type:

dict

riemannian.test_distance_and_exp_map()

Exponential map and geodesic distance (2D).

riemannian.test_error_handling_1d()

Invalid inputs raise appropriate errors (1D).

riemannian.test_euclidean_metric()

Flat Euclidean metric (2D).

riemannian.test_geodesic_integrators_2d()

Symplectic energy conservation and arc-length reparametrisation (2D).

riemannian.test_geodesic_integrators_accuracy()

Accuracy of 1D geodesic integrators on flat metric.

riemannian.test_hamiltonian_construction()

Metric from Hamiltonian (2D).

riemannian.test_hamiltonian_flow_conservation()

Energy conservation in 1D Hamiltonian flow.

riemannian.test_hodge_star_and_volume()

Hodge star and Riemannian volume (2D).

riemannian.test_integration_methods_consistency()

Symbolic vs numerical volume integration (1D).

riemannian.test_jacobi_field_stability()

Jacobi equation on the sphere (2D).

riemannian.test_laplace_beltrami_properties()

Laplace-Beltrami symbol (1D).

riemannian.test_metric_geometry_and_sl()

Geometric quantities and Sturm-Liouville reduction (1D).

riemannian.test_poincare_half_plane()

Poincaré half-plane — constant curvature -1 (2D).

riemannian.test_polar_coordinates()

Polar coordinate metric (2D).

riemannian.test_sphere_metric()

Unit sphere metric (2D).

riemannian.verify_gauss_bonnet(metric, domain, resolution=100)

Numerically verify the Gauss-Bonnet theorem (2D only).

∫∫_M K dA ≈ 2π χ(M)

Parameters:

• metric (Metric (dim=2))

• domain (tuple ((x_min, x_max), (y_min, y_max)))

• resolution (int)

Returns:

‘integral’, ‘expected’ (2π), ‘integration_error’, ‘relative_error’.

Return type:

dict

riemannian.visualize_curvature(metric, x_range=None, y_range=None, resolution=100, quantity='gauss', cmap='RdBu_r', **kwargs)

Visualize curvature of a 1D or 2D manifold.

Parameters:

• metric (Metric)

• x_range (tuple, optional)

• y_range (tuple, optional (2D only))

• resolution (int)

• quantity (str) – 2D: ‘gauss’ | ‘ricci_scalar’. 1D: ‘metric’ | ‘christoffel’.

• cmap (str)

• **kwargs (extra options for 1D (initial_conditions, tspan, colorby, n_steps).)

riemannian.visualize_geodesics(metric, initial_conditions, tspan, x_range=None, y_range=None, colorby='speed', plot_curvature=True, n_steps=500)

Visualize geodesics on a 1D or 2D manifold.

Parameters:

• metric (Metric)

• initial_conditions (list) – 1D: list of (x₀, v₀). 2D: list of ((x₀, y₀), (vₓ₀, vᵧ₀)).

• tspan (tuple)

• x_range (tuple, optional) – Plotting range (inferred if None).

• y_range (tuple, optional) – Plotting range (inferred if None).

• colorby ({'speed', 'time', 'curvature'}) – 1D only.

• plot_curvature (bool) – 2D only — show Gaussian curvature as background.

• n_steps (int)