riemannian

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

Overview

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

The module is organised in four broad layers:

1. Symbolic core — Metric and its methods build and cache all derived geometric quantities (inverse metric, Christoffel symbols, curvature tensors, Laplace–Beltrami symbol) as SymPy expressions.

2. Geodesic layer — standalone functions integrate geodesic and Hamiltonian flow, compute geodesic distance, solve the Jacobi deviation equation, and perform parallel transport along curves.

3. Differential-form layer (2D only) — the Hodge star, the de Rham Laplacian with Weitzenböck correction, and a full numerical Hodge decomposition backed by a sparse FEM solver on a RiemannianGrid.

4. Visualisation layer — geodesic trajectory plots and curvature colour maps for both dimensions.

Public API summary

Metric class

Central object; accepts a 1D scalar expression or a 2×2 SymPy Matrix together with coordinate symbols. Can also be constructed from a free-particle Hamiltonian via Metric.from_hamiltonian.

Methods available on every Metric:

• eval(*point) — evaluate all metric quantities at a numerical point.

• gauss_curvature() — Gaussian curvature K (2D) or zero (1D).

• riemann_tensor() — full Riemann curvature tensor R^i_jkl (2D).

• ricci_tensor() — Ricci tensor R_ij (2D).

• ricci_scalar() — scalar curvature R (2D).

• laplace_beltrami_symbol() — principal and subprincipal symbols of the Laplace-Beltrami operator, ready for microlocal analysis.

• riemannian_volume(domain) — Riemannian volume (area in 2D) by symbolic or numerical integration.

• arc_length(x_min, x_max) — arc length along a 1D metric.

• covariant_derivative_vector(V) — covariant derivative of a contravariant vector field (2D).

• covariant_derivative_covector(omega) — covariant derivative of a 1-form (2D).

• riemannian_gradient(f) — gradient of a scalar as a contravariant vector field, grad f = g^ij d_j f d_i (2D).

• riemannian_hessian(f) — covariant Hessian Hess(f)_ij = d_i d_j f - Gamma^k_ij d_k f (2D).

Geodesic and flow functions

• geodesic_solver(metric, p0, v0, tspan) — integrate the geodesic ODE with selectable schemes: 'rk4', 'adaptive', 'symplectic'.

• geodesic_hamiltonian_flow(metric, p0, v0, tspan) — geodesic flow in Hamiltonian form via the companion symplectic module (Verlet and other symplectic integrators).

• distance(metric, p, q) — geodesic distance by shooting or optimisation (2D).

• jacobi_equation_solver(metric, geodesic, initial_variation, tspan) — solve the Jacobi (geodesic deviation) equation along a reference geodesic (2D).

• parallel_transport(metric, curve, initial_vector, tspan) — parallel transport of a vector along a parameterised curve (2D).

• christoffel(metric) — convenience wrapper returning the Christoffel symbol dict.

• exponential_map(metric, p, v, t) — exponential map exp_p(tv) at a base point (2D).

Spectral and operator functions

• laplace_beltrami(metric) — return the Laplace-Beltrami symbol dict (wrapper around Metric.laplace_beltrami_symbol).

• sturm_liouville_reduce(metric, potential_expr) — reduce the 1D Laplace-Beltrami eigenvalue problem -Delta_g u + Vu = lambda u to canonical Sturm-Liouville form (1D only).

Differential-form layer (2D only)

• hodge_star(metric, form_degree) — symbolic Hodge star on 0-, 1-, and 2-forms.

• de_rham_laplacian(metric, form_degree) — symbolic de Rham-Hodge Laplacian Delta = d delta + delta d on 0- and 1-forms, including the Weitzenbock curvature correction K*id for 1-forms; returns a dict with the principal symbol, subprincipal symbol, Weitzenbock term, and a callable 'action' that applies the operator symbolically.

• hodge_decomposition(metric, omega, domain, form_degree) — numerical Hodge-Helmholtz decomposition of a 1-form or 2-form into exact, co-exact, and harmonic parts on a rectangular domain. The exact potential phi is solved with homogeneous Dirichlet BC (phi = 0 on dOmega) and the co-exact potential psi with Neumann-pinned BC (d_n psi = 0 on dOmega, gauge fixed at one interior node). The harmonic residual captures the de Rham cohomology class of the input form.

• visualize_hodge_decomposition(decomp) — vector-field or scalar-field plot of the three Hodge components with a metric-weighted energy bar chart.

RiemannianGrid class

Helper assembled once from a Metric and a rectangular domain. Caches the N^2 x N^2 sparse FEM matrix for Delta_0 (scalar Laplace-Beltrami) and the 2N^2 x 2N^2 block matrix for Delta_1 (de Rham Laplacian on 1-forms with Weitzenbock correction K*id). Exposes two public solvers:

• solve_poisson_dirichlet(rhs) — Delta_0 u = f with u = 0 on dOmega.

• solve_poisson_neumann(rhs) — Delta_0 u = f with d_n u = 0 on dOmega; the RHS is automatically projected onto the compatible subspace (mean-zero with respect to sqrt(g)), and the gauge is fixed by pinning one interior node to zero. Both solvers accept a single (N, N) array or a batch (K, N, N) array of right-hand sides.

Gauss-Bonnet verification

• verify_gauss_bonnet(metric, domain) — numerically integrate K dA over a rectangular domain and compare with the topological prediction 2 pi chi (2D only).

Numerical curvature helpers

• _brioschi_curvature_grid(g11, g12, g22, du, dv) — compute Gaussian curvature K on a meshgrid via the Brioschi formula, using only the metric components and their finite-difference derivatives. Used internally by build_embedding in place of the symbolic gauss_curvature() path, and available as a standalone helper for fast numerical K evaluation in the visualisation and corrugation layers.

Embedding and Nash–Kuiper corrugation layer** (2D only)

• build_embedding(metric, u_range, v_range, nu, nv) — constructs an approximate C¹ isometric embedding of a 2D metric into ℝ³ by marching row-by-row along the v-direction while maintaining an orthonormal Darboux frame (E1, E2, N) at each grid point. Gaussian curvature is computed numerically via _brioschi_curvature_grid and used to evolve the surface normal through the Weingarten map. The result is a first-order explicit approximation; the quality of the embedding can be assessed with metric_deficit.

• metric_deficit(R, g11_t, g12_t, g22_t, du, dv) — computes the pointwise difference between the target metric components and the metric induced by the embedding R, returning the three deficit grids (Δg₁₁, Δg₁₂, Δg₂₂) together with the overall Frobenius norm of the deficit.

• corrugation_step(R, dg11, dg12, dg22, du, dv, freq) — applies one Nash–Kuiper corrugation pass: the deficit tensor is decomposed into rank-1 eigenvalue terms at each grid point, and sinusoidal normal oscillations of amplitude ρ = √λ / (√2 π freq) are added so that each positive eigenvalue is filled in the phase-averaged sense.

• add_corrugations(R, metric, du, dv, U, V, n_iterations, base_freq, alpha) — orchestrates the full Nash–Kuiper pipeline: it first scales R to a short map with factor α < 1, then repeatedly calls corrugation_step with geometrically doubling spatial frequency (base_freq · 2ⁱ) to drive the Frobenius deficit toward zero. Returns a dict with the short map, the final embedding, per-iteration snapshots, and the deficit convergence history.

• plot_embedding and plot_corrugation_pipeline complete the embedding workflow with 3D surface visualisation and a convergence diagnostic figure, respectively.

Visualisation

• visualize_geodesics(metric, initial_conditions, tspan) — overlay geodesic trajectories on a curvature background (1D or 2D).

• visualize_curvature(metric, x_range, y_range, quantity) — colour map of Gaussian curvature or Ricci scalar (2D), or line plot of the metric component / Christoffel symbol with optional geodesic overlay (1D).

Mathematical background

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

ds^2 = g_ij(x) dx^i dx^j,

with inverse g^ij. The Christoffel symbols are

Gamma^i_jk = 1/2 * g^il (d_j g_kl + d_k g_jl - d_l g_jk)

and govern the geodesic equation x^i_tt + Gamma^i_jk x^j_t x^k_t = 0. In 1D this simplifies to x_tt + Gamma^1_11 x_t^2 = 0 with Gamma^1_11 = 1/2 (log g_11)’.

Curvature is encoded in the Riemann tensor R^i_jkl, the Ricci tensor R_ij = R^k_ikj, and the scalar curvature R = g^ij R_ij. In 2D the Gaussian curvature K is the only independent component: R_1212 = K |g|.

The Laplace-Beltrami operator on scalar functions is

Delta = |g|^{-1/2} d_i ( |g|^{1/2} g^ij d_j ),

with principal symbol g^ij xi_i xi_j. The sparse FEM discretisation uses a standard 5-point finite-volume stencil with cross-derivative terms for off-diagonal metric components, divided by sqrt(g) to recover the true operator.

The de Rham Laplacian Delta = d delta + delta d on k-forms coincides with the Laplace-Beltrami operator on 0-forms and satisfies the Weitzenbock identity on 1-forms in 2D:

Delta alpha = nabla* nabla alpha + K * alpha,

where nabla* nabla is the rough (connection) Laplacian acting component-wise and K is the Gaussian curvature.

The Hodge decomposition of a differential form on a compact domain Omega in 2D is as follows:

• 0‑form (scalar function) f:

  f = Δφ + h₀, where Δφ = δ(⋆dφ) is the co‑exact part (the Laplacian of a scalar potential φ) and h₀ is the constant harmonic part (the mean of f). The exact part is zero because no (−1)-form exists. The harmonic part corresponds to the de Rham cohomology class [f] ∈ H⁰_dR(Ω) ≅ ℝ.

• 1‑form α = αₓ dx + αᵧ dy:

  α = dφ + ⋆dψ + h, where dφ is exact, ⋆dψ is co‑exact, and h is harmonic (Δh = 0). φ and ψ are scalar potentials satisfying Poisson equations; h represents the de Rham cohomology class [α] ∈ H¹_dR(Ω).

• 2‑form ω = f dx∧dy:

  ω = d(⋆dφ) + h₂, with no co‑exact component (no 3‑forms exist in 2D). φ is a scalar potential solving Δφ = -⋆ω, and h₂ is the harmonic 2‑form, representing the cohomology class [ω] ∈ H²_dR(Ω).

References

[1]

do Carmo, M. P. Riemannian Geometry, Birkhauser, 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.).

[6]

Warner, F. W. Foundations of Differentiable Manifolds and Lie Groups, Springer, 1983.

class riemannian.Metric(g_input, coords)

Bases: object

Riemannian metric on a 1D or 2D manifold.

The dimension is inferred automatically from the supplied inputs:

• 1D: g_input is a scalar SymPy expression in one coordinate; coords is a 1-tuple (x,) or a bare symbol.

• 2D: g_input is a 2×2 SymPy Matrix (or a nested list that will be promoted to one); coords is a 2-tuple (x, y).

On construction the metric is simplified once, and all derived symbolic quantities — inverse metric, determinant, square-root of the determinant, and Christoffel symbols — are computed and stored. Numerical callables (lambdify-produced functions) are also built and cached so that integration and visualization routines have zero symbolic overhead at run-time.

Parameters:

• g_input (sympy.Expr or sympy.Matrix or list) – Metric tensor. A scalar SymPy expression for 1D, or a 2×2 SymPy Matrix (or equivalent nested list of expressions) for 2D.

• coords (tuple of sympy.Symbol or sympy.Symbol) – Coordinate symbols in order. The length determines the manifold dimension (1 or 2). A bare symbol is accepted for the 1D case.

dim

Manifold dimension, either 1 or 2.

Type:

int

coords

Coordinate symbols, always stored as a tuple.

Type:

tuple of sympy.Symbol

g_expr

(1D only) Simplified metric component g₁₁(x).

Type:

sympy.Expr

g_inv_expr

(1D only) Inverse metric component g¹¹(x) = 1/g₁₁(x).

Type:

sympy.Expr

sqrt_det_expr

(1D only) Square root of the metric determinant, √|g₁₁(x)|.

Type:

sympy.Expr

christoffel_sym

Symbolic Christoffel symbols. - 1D: Expression for Γ¹₁₁ = ½ (log g₁₁)′. - 2D: Nested dict christoffel_sym[i][j][k] → SymPy expression.

Type:

sympy.Expr or dict

g_func

Numerical metric function. - 1D: g₁₁(x_val) → float or ndarray. - 2D: Dict {(i, j): callable} of component functions.

Type:

callable or dict

g_inv_func

Numerical inverse-metric function. - 1D: g¹¹(x_val) → float or ndarray. - 2D: Dict {(i, j): callable} of component functions.

Type:

callable or dict

sqrt_det_func

(1D only) Numerical function √|g₁₁|(x_val) → float or ndarray.

Type:

callable

sqrt_det_g_func

(2D only) Numerical function √|det(g)|(x_val, y_val).

Type:

callable

christoffel_func

Numerical Christoffel callables. - 1D: Γ¹₁₁(x_val) → float or ndarray. - 2D: Nested dict of callables Γⁱⱼₖ(x_val, y_val).

Type:

callable or dict

g_matrix

(2D only) Simplified 2×2 metric tensor matrix.

Type:

sympy.Matrix

det_g

(2D only) Determinant of the metric, det(g).

Type:

sympy.Expr

sqrt_det_g

(2D only) Square root of the absolute determinant, √|det(g)|.

Type:

sympy.Expr

g_inv_matrix

(2D only) Symbolic inverse metric g⁻¹.

Type:

sympy.Matrix

det_g_func

(2D only) Numerical function det(g)(x_val, y_val).

Type:

callable

Raises:

• ValueError – If the number of coordinate symbols is neither 1 nor 2.

• ValueError – If g_input is not a 2×2 matrix when len(coords) == 2.

Examples

1D — cone-like metric g = x²:

>>> from sympy import symbols, Matrix, sin, simplify
>>> x = symbols('x', real=True, positive=True)
>>> m = Metric(x**2, (x,))
>>> m.dim
1
>>> m.christoffel_sym          # Γ¹₁₁ = 1/x
1/x
>>> m.gauss_curvature()        # intrinsic curvature of a curve is 0
0

2D — unit sphere ds² = dθ² + sin²θ dφ²:

>>> theta, phi = symbols('theta phi', real=True)
>>> g = Matrix([[1, 0], [0, sin(theta)**2]])
>>> m = Metric(g, (theta, phi))
>>> m.dim
2
>>> simplify(m.gauss_curvature())   # K = 1 everywhere on the unit sphere
1

2D — Poincaré half-plane g = diag(1/y², 1/y²):

>>> x, y = symbols('x y', real=True)
>>> m = Metric(Matrix([[1/y**2, 0], [0, 1/y**2]]), (x, y))
>>> simplify(m.gauss_curvature())   # constant negative curvature
-1

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

Compute the arc length of a 1D metric between two coordinate values.

Convenience wrapper around riemannian_volume() for the 1D case. The arc length is

L = ∫_{x_min}^{x_max} √g₁₁(x) dx.

Parameters:

• x_min (float or sympy.Expr) – Integration bounds along the coordinate x.

• x_max (float or sympy.Expr) – Integration bounds along the coordinate x.

• method ({'symbolic', 'numerical'}, default 'numerical') – Passed directly to riemannian_volume().

Returns:

Arc length, symbolic or numerical depending on method.

Return type:

sympy.Expr or float

Raises:

NotImplementedError – If called on a 2D metric. Use riemannian_volume() with a 2D domain for surface area computations.

covariant_derivative_covector(omega_components, do_simplify=True)

Compute the covariant derivative ∇ω of a covector field (2D only).

The covariant derivative of a covariant 1-form ω = ωⱼ dxʲ in the direction of the coordinate basis vector ∂ᵢ is

(∇ω)ᵢⱼ = ∂ᵢ ωⱼ − Γᵏᵢⱼ ωₖ,

where the Christoffel term accounts for the change of basis under parallel transport.

Parameters:

• omega_components (list or tuple of two sympy.Expr) – Lower-index components (ω₁, ω₂) of the 1-form, expressed as functions of the coordinate symbols.

• do_simplify (bool, default True) – If True, each entry of the resulting matrix is passed through sympy.simplify before being returned. Set to False for faster (but potentially unsimplified) output.

Returns:

A 2×2 matrix whose (i, j)-entry is ∇ᵢ ωⱼ, i.e. the covariant derivative of the j-th component of ω in the i-th coordinate direction. Both row index i and column index j are lower (covariant) indices taking values in {0, 1}.

Return type:

sympy.Matrix

Raises:

NotImplementedError – If called on a 1D metric.

Examples

>>> from sympy import symbols, Matrix
>>> x, y = symbols('x y', real=True)
>>> m = Metric(Matrix([[1, 0], [0, 1]]), (x, y))   # flat metric
>>> nabla_omega = m.covariant_derivative_covector([x**2, y**2])
>>> nabla_omega   # Christoffel terms vanish on flat space
Matrix([[2*x, 0], [0, 2*y]])

covariant_derivative_vector(V_components, do_simplify=True)

Compute the covariant derivative ∇V of a vector field (2D only).

The covariant derivative of a contravariant vector field V = Vʲ ∂ⱼ in the direction of the coordinate basis vector ∂ᵢ is

(∇V)ⁱⱼ = ∂ᵢ Vʲ + Γʲᵢₖ Vᵏ,

where Γʲᵢₖ are the Christoffel symbols of the second kind.

Parameters:

• V_components (list or tuple of two sympy.Expr) – Upper-index components (V¹, V²) of the vector field, expressed as functions of the coordinate symbols.

• do_simplify (bool, default True) – If True, each entry of the resulting matrix is passed through sympy.simplify before being returned. Set to False for faster (but potentially unsimplified) output.

Returns:

A 2×2 matrix whose (i, j)-entry is ∇ᵢ Vʲ, i.e. the j-th component of the covariant derivative in the i-th coordinate direction. Row index i = 0, 1 is the differentiation direction (lower); column index j = 0, 1 is the vector component (upper).

Return type:

sympy.Matrix

Raises:

NotImplementedError – If called on a 1D metric.

Examples

>>> from sympy import symbols, Matrix
>>> x, y = symbols('x y', real=True)
>>> m = Metric(Matrix([[1, 0], [0, 1]]), (x, y))   # flat metric
>>> nabla_V = m.covariant_derivative_vector([x, y])
>>> nabla_V   # equals the ordinary Jacobian for flat space
Matrix([[1, 0], [0, 1]])

eval(*point)

Evaluate metric quantities numerically at a given coordinate point.

All cached numerical callables are evaluated at the supplied point and collected into a single dict for convenient downstream use (e.g. inside ODE right-hand sides or distance computations).

Parameters:

*point (float or numpy.ndarray) – Coordinate values at which to evaluate. Supply one argument for a 1D metric and two arguments for a 2D metric. NumPy arrays are accepted, in which case all returned values are arrays of the same shape.

Returns:

1D — keys and values:

• 'g' : float or ndarray — metric component g₁₁(x).

• 'g_inv' : float or ndarray — inverse metric g¹¹(x).

• 'sqrt_det' : float or ndarray — √|g₁₁(x)|.

• 'christoffel' : float or ndarray — Γ¹₁₁(x).

2D — keys and values:

• 'g' : ndarray of shape (2, 2, …) — full metric matrix.

• 'g_inv' : ndarray of shape (2, 2, …) — inverse metric.

• 'det_g' : float or ndarray — det(g).

• 'sqrt_det' : float or ndarray — √|det(g)|.

• 'christoffel' : nested dict [i][j][k] — Γⁱⱼₖ values.

Return type:

dict

Examples

>>> x = symbols('x', real=True, positive=True)
>>> m = Metric(x**2, (x,))
>>> ev = m.eval(2.0)
>>> ev['g']             # g₁₁(2) = 4
4.0
>>> ev['christoffel']   # Γ¹₁₁(2) = 1/2
0.5

classmethod from_hamiltonian(H_expr, coords, momenta)

Construct a Metric by extracting the kinetic term from a Hamiltonian.

For a free-particle (purely kinetic) Hamiltonian of the form

H = ½ gⁱʲ(q) pᵢ pⱼ + V(q),

the contravariant metric tensor gⁱʲ equals the Hessian of H with respect to the momenta. This method computes that Hessian and inverts it to recover the covariant metric gᵢⱼ. Any potential term V(q) that is independent of the momenta is automatically discarded.

Parameters:

• H_expr (sympy.Expr) – Full Hamiltonian expressed in terms of both coords and momenta symbols.

• coords (tuple of sympy.Symbol) – Generalised position variables, length 1 (1D) or 2 (2D).

• momenta (tuple of sympy.Symbol) – Conjugate momentum variables, same length as coords.

Returns:

A Metric instance whose dimension matches len(coords).

Return type:

Metric

Raises:

ValueError – If len(coords) is neither 1 nor 2.

Notes

The extraction relies on ∂²H/∂pᵢ∂pⱼ = gⁱʲ, which is exact for quadratic kinetic terms. Non-quadratic kinetic energy (e.g. relativistic Hamiltonians) will produce incorrect results.

Examples

1D polar-like kinetic energy H = p² / (2x²):

>>> from sympy import symbols, simplify
>>> x, p = symbols('x p', real=True)
>>> H = p**2 / (2*x**2)
>>> m = Metric.from_hamiltonian(H, (x,), (p,))
>>> m.dim
1
>>> simplify(m.g_expr - x**2)    # recovers g₁₁ = x²
0

2D polar coordinates H = (p_r² + p_θ²/r²) / 2:

>>> r, theta = symbols('r theta', real=True, positive=True)
>>> pr, pt = symbols('p_r p_theta', real=True)
>>> H = (pr**2 + pt**2/r**2) / 2
>>> m = Metric.from_hamiltonian(H, (r, theta), (pr, pt))
>>> m.dim
2

gauss_curvature()

Compute the Gaussian curvature K of the manifold.

For a 1D manifold (a curve), the intrinsic Gaussian curvature is identically zero and the method returns the SymPy integer 0 immediately without any computation.

For a 2D surface the curvature is extracted from the Riemann tensor via

K = R₁₂₁₂ / det(g),

where R₁₂₁₂ = gₐ₁ Rᵃ₂₁₂ is the single independent component of the fully covariant Riemann tensor in 2D.

Returns:

Simplified symbolic expression for K as a function of the coordinates. For flat metrics this simplifies to 0; for the unit sphere it simplifies to 1.

Return type:

sympy.Expr

Notes

The Riemann tensor is recomputed each call (no caching). For repeated curvature queries on the same metric consider caching the result yourself, e.g. K = simplify(m.gauss_curvature()).

Examples

>>> from sympy import symbols, Matrix, sin, simplify
>>> x, y = symbols('x y', real=True)
>>> m_flat = Metric(Matrix([[1, 0], [0, 1]]), (x, y))
>>> m_flat.gauss_curvature()
0

>>> theta, phi = symbols('theta phi', real=True)
>>> m_sphere = Metric(Matrix([[1, 0], [0, sin(theta)**2]]), (theta, phi))
>>> simplify(m_sphere.gauss_curvature())
1

laplace_beltrami_symbol()

Compute the full symbol of the Laplace–Beltrami operator.

The Laplace–Beltrami operator acting on smooth functions is

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

Its principal symbol (the highest-order part) evaluated at cotangent vector ξ is the positive-definite quadratic form

σ₂(Δ_g)(x, ξ) = gⁱʲ(x) ξᵢ ξⱼ.

The subprincipal symbol captures the first-order (transport) part arising from the non-constancy of the metric. The two are combined into a full symbol via

σ_full = σ₂ + i σ₁ (microlocal convention),

making the output ready for use in microlocal analysis or WKB approximations.

Works for both 1D and 2D metrics.

Returns:

• dict with keys

• * ``’principal’`` (sympy.Expr — σ₂(Δ_g)(x, ξ), a polynomial of) – degree 2 in the cotangent variables ξ (and η for 2D).

• * ``’subprincipal’`` (sympy.Expr — σ₁(Δ_g)(x, ξ), degree 1 in ξ.)

• * ``’full’`` (sympy.Expr — σ₂ + i σ₁ (complex-valued symbol).)

Notes

The cotangent variables are introduced as fresh SymPy symbols named xi (1D) or xi, eta (2D). They do not conflict with the coordinate symbols.

Examples

1D, g = x²:

>>> from sympy import symbols, simplify
>>> x, xi = symbols('x xi', real=True, positive=True)
>>> m = Metric(x**2, (x,))
>>> lb = m.laplace_beltrami_symbol()
>>> simplify(lb['principal'] - xi**2/x**2)
0
>>> simplify(lb['subprincipal'] - xi/x**3)
0

2D, polar coordinates g = diag(1, r²):

>>> r, theta = symbols('r theta', real=True, positive=True)
>>> xi_r, xi_t = symbols('xi eta', real=True)
>>> m = Metric(Matrix([[1, 0], [0, r**2]]), (r, theta))
>>> lb = m.laplace_beltrami_symbol()
>>> simplify(lb['principal'] - (xi_r**2 + xi_t**2/r**2))
0

ricci_scalar()

Compute the scalar (Ricci) curvature R.

The scalar curvature is the full trace of the Ricci tensor with respect to the metric:

R = gⁱʲ Rᵢⱼ.

For 1D manifolds the scalar curvature is identically zero and the method returns the SymPy integer 0 immediately. For 2D surfaces R = 2K, where K is the Gaussian curvature.

Returns:

Simplified symbolic expression for the scalar curvature.

Return type:

sympy.Expr

Notes

This method internally calls ricci_tensor(), which in turn calls riemann_tensor(). Both are recomputed on each invocation.

Examples

>>> from sympy import symbols, Matrix, sin, simplify
>>> theta, phi = symbols('theta phi', real=True)
>>> m = Metric(Matrix([[1, 0], [0, sin(theta)**2]]), (theta, phi))
>>> simplify(m.ricci_scalar())   # R = 2 for the unit sphere
2

ricci_tensor()

Compute the Ricci tensor Rᵢⱼ (2D manifolds only).

The Ricci tensor is the contraction

Rᵢⱼ = Rᵏᵢₖⱼ

of the Riemann tensor over the first and third indices. In 2D the result is a symmetric 2×2 SymPy matrix.

Returns:

Simplified 2×2 symmetric matrix of Ricci tensor components.

Return type:

sympy.Matrix

Raises:

NotImplementedError – If called on a 1D metric.

Examples

>>> from sympy import symbols, Matrix, sin, simplify
>>> theta, phi = symbols('theta phi', real=True)
>>> m = Metric(Matrix([[1, 0], [0, sin(theta)**2]]), (theta, phi))
>>> Ric = m.ricci_tensor()
>>> simplify(Ric[0, 0])   # R_θθ = 1 on the unit sphere
1

riemann_tensor()

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

The tensor is computed from the Christoffel symbols via the standard formula

Rⁱⱼₖₗ = ∂ₖΓⁱⱼₗ − ∂ₗΓⁱⱼₖ + ΓⁱₘₖΓᵐⱼₗ − ΓⁱₘₗΓᵐⱼₖ,

with all 16 components evaluated and simplified. In 2D the tensor has at most one independent component (R⁰₁₀₁ = −R⁰₁₁₀ = R¹₀₁₀ = …), and the Gaussian curvature is proportional to it.

Returns:

Nested dict R[i][j][k][l] → SymPy expression for Rⁱⱼₖₗ, with all indices in {0, 1}.

Return type:

dict

Raises:

NotImplementedError – If called on a 1D metric (the Riemann tensor of a 1D manifold is identically zero and carries no information).

Notes

All 16 components are computed and returned; the user is responsible for exploiting symmetries (antisymmetry in k, l and in i, j when lowered) to reduce redundant evaluations.

Examples

>>> from sympy import symbols, Matrix, sin, simplify
>>> theta, phi = symbols('theta phi', real=True)
>>> m = Metric(Matrix([[1, 0], [0, sin(theta)**2]]), (theta, phi))
>>> R = m.riemann_tensor()
>>> simplify(R[0][1][0][1])   # non-zero component on unit sphere
-sin(theta)**2

riemannian_gradient(f_expr, do_simplify=True)

Compute the gradient of a scalar function f as a contravariant vector field.

In local coordinates: ∇f = gⁱʲ ∂ⱼ f ∂ᵢ.

Parameters:

• f_expr (sympy.Expr) – The scalar function expressed in the metric’s coordinates.

• do_simplify (bool, default True) – Whether to simplify the resulting expressions.

Returns:

• 1D: a single SymPy expression for the gradient component (∇f)¹.

• 2D: a tuple (∇f¹, ∇f²) of SymPy expressions.

Return type:

tuple or dict

riemannian_hessian(f_expr, do_simplify=True)

Compute the Hessian of a scalar function f as a (0,2)-tensor (covariant).

In local coordinates: Hess(f)ᵢⱼ = ∂ᵢ∂ⱼ f − Γᵏᵢⱼ ∂ₖ f.

Parameters:

• f_expr (sympy.Expr) – The scalar function expressed in the metric’s coordinates.

• do_simplify (bool, default True) – Whether to simplify the resulting expressions.

Returns:

The Hessian matrix (2×2) for 2D, or the single component H₁₁ for 1D.

Return type:

sympy.Matrix (2D) or sympy.Expr (1D)

riemannian_volume(domain, method='numerical')

Compute the Riemannian volume (area / arc length) of a rectangular domain.

The Riemannian volume element is dV = √|det(g)| dx¹ ∧ … ∧ dxⁿ. For a 1D metric this equals the arc length ∫ √g₁₁ dx; for a 2D metric it equals the surface area ∫∫ √|det(g)| dx dy.

Parameters:

• domain (tuple) –

  • 1D: (x_min, x_max) — a pair of floats or SymPy numbers defining the integration interval.

  • 2D: ((x_min, x_max), (y_min, y_max)) — a pair of pairs defining the rectangular integration region.

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

  • 'symbolic': uses SymPy’s integrate for an exact result. May be slow or fail to close for complicated metrics.

  • 'numerical': uses scipy.integrate.quad (1D) or scipy.integrate.dblquad (2D) for a fast floating-point approximation.

Returns:

• Symbolic result (SymPy expression or number) when method='symbolic'.

• Float when method='numerical'.

Return type:

sympy.Expr or float

Raises:

ValueError – If method is not 'symbolic' or 'numerical'.

Examples

1D — arc length of g = 1/x² on [1, e]:

>>> import numpy as np
>>> x = symbols('x', real=True, positive=True)
>>> m = Metric(1/x**2, (x,))
>>> float(m.riemannian_volume((1, np.e), method='symbolic'))  # = 1
1.0
>>> m.riemannian_volume((1, np.e), method='numerical')         # ≈ 1
0.9999...

2D — area of g = diag(4, 9) on [0,1]²:

>>> x, y = symbols('x y', real=True)
>>> m = Metric(Matrix([[4, 0], [0, 9]]), (x, y))
>>> m.riemannian_volume(((0, 1), (0, 1)), method='symbolic')   # = 6
6

class riemannian.RiemannianGrid(metric, domain, resolution)

Bases: object

A regular rectangular grid carrying numerical metric data and the assembled sparse FEM Laplace–Beltrami matrix.

This is the bridge between the symbolic de_rham_laplacian operator and the numerical solvers in hodge_decomposition. Instead of each function building its own stencil, both consume a RiemannianGrid.

Parameters:

• metric (Metric) – A 2D Riemannian metric.

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

• resolution (int) – Number of grid points per axis.

X, Y

Meshgrid coordinate arrays.

Type:

ndarray (N, N)

dx, dy

Grid spacings.

Type:

float

N, N2

Points per axis and total points.

Type:

int

sqrt_det, g_inv00, g_inv11, g_inv01

Evaluated metric coefficients.

Type:

ndarray (N, N)

A_scalar

Assembled FEM matrix for the scalar Laplace–Beltrami operator Δ₀. Suitable for solving Δ₀ u = f after applying Dirichlet BC.

Type:

scipy.sparse.csr_matrix

A_1form

Block-2×2 matrix for the de Rham Laplacian on 1-forms, incorporating the Weitzenböck curvature term. Shape (2N², 2N²).

Type:

scipy.sparse.csr_matrix

idx_bound

Flat indices of all boundary nodes (for Dirichlet BC enforcement).

Type:

ndarray of int

Examples

>>> from sympy import symbols, Matrix
>>> x, y = symbols('x y', real=True)
>>> m = Metric(Matrix([[1, 0], [0, 1]]), (x, y))
>>> grid = RiemannianGrid(m, ((0, 1), (0, 1)), resolution=30)
>>> grid.A_scalar.shape
(900, 900)
>>> grid.A_1form.shape
(1800, 1800)

solve_poisson_dirichlet(rhs)

Solve the scalar Laplace–Beltrami equation Δ₀ u = rhs with homogeneous Dirichlet boundary conditions.

The system is assembled as a sparse FEM matrix self._A_dirichlet, which encodes the operator Δ₀ on a regular rectangular grid. Homogeneous Dirichlet conditions (u = 0 on the boundary) are enforced by zeroing the rows corresponding to boundary nodes and setting the diagonal entry to 1; the right‑hand side is set to zero on those rows.

The solver supports multiple right‑hand sides simultaneously, which is more efficient than solving each individually.

Parameters:

rhs (numpy.ndarray) –

Right‑hand side field(s). Allowed shapes:

• 2D (N, N) : a single scalar field on the grid.

• 3D (K, N, N) : K independent scalar fields, solved in one linear system.

Returns:

Solution u of the Poisson equation, with the same shape as rhs:

• (N, N) for a single RHS.

• (K, N, N) for multiple RHS (first dimension = number of fields).

Return type:

numpy.ndarray

Raises:

ValueError – If the shape of rhs does not match the grid dimensions.

Notes

The linear system is solved using scipy.sparse.linalg.spsolve, which returns the solution with the same number of columns as the number of RHS vectors. The result is reshaped back to the original spatial grid.

solve_poisson_neumann(rhs)

Solve the scalar Laplace–Beltrami equation Δ₀ u = rhs with homogeneous Neumann boundary conditions.

The Neumann problem requires a compatibility condition: the right‑hand side must be orthogonal to the constant nullspace, i.e. ∫_M rhs · √g dV = 0. This function automatically enforces that condition by subtracting the weighted mean from rhs before solving. The gauge freedom (addition of an arbitrary constant) is fixed by pinning a single interior node (the grid point nearest to the centre) to zero in the assembled matrix self._A_neumann.

The solver supports multiple right‑hand sides simultaneously, which is more efficient than solving each individually.

Parameters:

rhs (numpy.ndarray) –

Right‑hand side field(s). Allowed shapes:

• 2D (N, N) : a single scalar field on the grid.

• 3D (K, N, N) : K independent scalar fields, solved in one linear system.

Returns:

Solution u of the Poisson equation, with the same shape as rhs:

• (N, N) for a single RHS.

• (K, N, N) for multiple RHS (first dimension = number of fields).

Return type:

numpy.ndarray

Raises:

ValueError – If the shape of rhs does not match the grid dimensions.

Notes

The weighted mean is computed as (∫ rhs · √g dV) / (∫ √g dV) and subtracted from each RHS vector. The solution is unique because the gauge is fixed by the interior node condition.

riemannian.add_corrugations(R, metric, du, dv, U, V, n_iterations=6, base_freq=2, alpha=0.85)

Iteratively apply Nash–Kuiper corrugations to drive the metric deficit to zero.

Workflow:

1. Short map: Scale R by alpha < 1 so that the induced metric strictly undershoots the target everywhere (a necessary precondition for the Nash–Kuiper scheme).

2. Corrugation loop: For each iteration k (0-indexed), compute the current metric deficit, apply corrugation_step at frequency base_freq · 2^k, and record the new Frobenius deficit. Doubling the frequency at each step ensures that successive passes add wrinkles at finer and finer scales without interfering destructively with earlier corrections.

Progress is printed to stdout at each iteration.

Parameters:

• R (ndarray, shape (nu, nv, 3)) – Initial embedding, typically produced by build_embedding.

• metric (Metric) – The target 2D Metric instance.

• du (float) – Uniform grid spacing in the u-direction.

• dv (float) – Uniform grid spacing in the v-direction.

• U (ndarray, shape (nu, nv)) – Meshgrid of u-coordinate values (indexing='ij').

• V (ndarray, shape (nu, nv)) – Meshgrid of v-coordinate values.

• n_iterations (int, optional) – Number of corrugation passes. Each pass halves the length scale of the wrinkles; 6–10 iterations are typical for smooth metrics. Default is 6.

• base_freq (float, optional) – Starting spatial frequency for the first corrugation pass. Frequencies used are base_freq · 2^k for k = 0, …, n_iterations-1. Default is 2.

• alpha (float, optional) – Short-map scaling factor; must be strictly less than 1. Values close to 1 (e.g. 0.85–0.95) minimise the initial deficit while still guaranteeing a true short map. Default is 0.85.

Returns:

result –

'R_short'ndarray, shape (nu, nv, 3)

The scaled short map used as the starting point (R · alpha).

'R_final'ndarray, shape (nu, nv, 3)

The embedding after all corrugation passes.

'snapshots'list of ndarray

One copy of R per iteration, starting from the short map (snapshots[0]) through the final result (snapshots[n_iterations]). Useful for animation.

'deficits'list of float

Frobenius norm of the metric deficit at each snapshot, in the same order as snapshots.

Return type:

dict with the following keys:

riemannian.analyze_hodge_decomposition(decomp, original=None, print_report=True, show_plot=True)

Analyse the output of hodge_decomposition and print key metrics.

Parameters:

• decomp (dict) – Dictionary returned by hodge_decomposition (must contain a ‘grid’ key).

• original (optional) – The original differential form used as input to hodge_decomposition. - For 1‑form: tuple of two arrays, callables, or SymPy expressions. - For 2‑form: a single array, callable, or SymPy expression. If not provided, reconstruction errors are not computed.

• print_report (bool, default True) – If True, print the analysis to the console.

• show_plot (bool, default True) – If True, call visualize_hodge_decomposition on decomp.

Returns:

Dictionary containing all computed metrics.

Return type:

dict

riemannian.build_embedding(metric, u_range, v_range, nu, nv)

Construct a C¹ 3D embedding that approximately realises a given 2D metric.

The algorithm marches row by row in the v-direction, maintaining an orthonormal Darboux frame (E1, E2, N) at each grid point:

• Initialisation (j = 0): The bottom row is placed along the x-axis with arc-length spacing determined by √g₁₁ · du. E1 points along x, N along z, and E2 along y. The v-derivative dR/dv is initialised from the metric components g₁₂ and g₂₂.

• Row advance (j → j+1):

  1. Positions are stepped as R[:, j+1] = R[:, j] + dv · dR/dv[:, j].

  2. The surface normal N is evolved via the Weingarten map using the shape-operator coefficient h₂₂, derived from K and the metric.

  3. E1 is re-estimated from a finite difference of the new row.

  4. E2 is recomputed as N × E1 to maintain right-handedness.

  5. dR/dv for the new row is reconstructed from g₁₂, g₂₂, E1, E2.

This is a first-order explicit integration scheme; accuracy improves with finer grids (larger nu, nv).

Parameters:

• metric (Metric) – A 2D Metric instance providing g_func. Gaussian curvature is now derived numerically via the Brioschi formula applied to the discrete metric grid, so gauss_curvature() is no longer called during embedding construction.

• u_range (tuple of float) – (u_min, u_max) — parameter range in the u-direction.

• v_range (tuple of float) – (v_min, v_max) — parameter range in the v-direction.

• nu (int) – Number of grid points in the u-direction.

• nv (int) – Number of grid points in the v-direction.

Returns:

• R (ndarray, shape (nu, nv, 3)) – The constructed embedding.

• u_vals (ndarray, shape (nu,)) – Uniformly spaced u parameter values.

• v_vals (ndarray, shape (nv,)) – Uniformly spaced v parameter values.

Notes

The embedding is not guaranteed to be isometric — use metric_deficit to assess how well g_induced matches g_target, and pass the result to add_corrugations to reduce the deficit iteratively.

The sign of the Gaussian curvature K determines whether h₂₂ causes the normal to tilt toward or away from the surface, which is how positive and negative curvature are distinguished geometrically.

riemannian.christoffel(metric)

Return the pre-computed numerical Christoffel symbol callable(s) of a metric.

This is a thin convenience accessor that exposes metric.christoffel_func without needing to remember the attribute name.

Parameters:

metric (Metric) – A Metric instance of any dimension.

Returns:

• 1D: a single callable Gamma(x_val) returning the scalar Γ¹₁₁(x) as a float or ndarray.

• 2D: a nested dict Gamma[i][j][k] of callables, each mapping (x_val, y_val) to Γⁱⱼₖ as a float or ndarray.

Return type:

callable or dict

Examples

>>> x = symbols('x', real=True, positive=True)
>>> m = Metric(x**2, (x,))
>>> G = christoffel(m)
>>> G(2.0)    # Γ¹₁₁(2) = 1/2
0.5

riemannian.corrugation_step(R, dg11, dg12, dg22, du, dv, freq)

Apply one Nash–Kuiper corrugation pass to reduce the metric deficit.

The Nash–Kuiper theorem guarantees that any short map into ℝ³ can be approximated arbitrarily closely by a C¹ isometric embedding via an infinite sequence of corrugation steps. Each step decomposes the symmetric deficit tensor Δg pointwise into two rank-1 terms via its eigendecomposition, then adds a sinusoidal normal oscillation that fills each positive eigenvalue in the phase-averaged sense.

For eigenvalue λ and eigenvector direction w, the oscillation is

R ← R + ρ · sin(2π·freq·⟨w, coord⟩ + φ) · n̂,

with amplitude ρ = √λ / (√2 · π · freq), chosen so that the added corrugation contributes exactly λ·w⊗w to the induced metric after phase averaging. Two phase offsets (0 and π/2) are applied per eigenvector to improve uniformity across the grid.

Parameters:

• R (ndarray, shape (N, N, 3)) – Current embedding. Must be square (nu == nv == N) as the coordinate phase arrays are built from a single size N.

• dg11 (ndarray, shape (N, N)) – Metric deficit in the g_{uu} component (target minus induced).

• dg12 (ndarray, shape (N, N)) – Metric deficit in the g_{uv} component.

• dg22 (ndarray, shape (N, N)) – Metric deficit in the g_{vv} component.

• du (float) – Grid spacing in the u-direction (used to scale the phase coordinate).

• dv (float) – Grid spacing in the v-direction.

• freq (float) – Spatial frequency of the corrugation oscillation. Higher values produce finer wrinkles. The corrugation formula is exact only in the limit freq → ∞; in practice freq should satisfy freq ≪ N / (2·L) (Nyquist limit relative to the domain size L). Recommended range: 2–8 for fine grids.

Returns:

R_new – Updated embedding after the corrugation pass. The input R is not modified in place.

Return type:

ndarray, shape (N, N, 3)

Notes

Only positive eigenvalues of Δg are corrugated; negative eigenvalues (where the induced metric already overshoots the target) are ignored. Repeatedly calling this function with doubling frequency (as done by add_corrugations) converges to an isometric embedding.

riemannian.de_rham_laplacian(metric, form_degree)

Compute the symbol and symbolic action of the de Rham–Hodge Laplacian Δ = dδ + δd on k-forms (2D only).

For form_degree=0 the result is identical to Metric.laplace_beltrami_symbol().

For form_degree=1 the Weitzenböck identity

Δα = ∇*∇α + Ric(α♯)♭

is used. In 2D the Ricci endomorphism on 1-forms is simply K·id (multiplication by the Gaussian curvature), so

(Δα)ᵢ = (∇*∇α)ᵢ + K αᵢ ,

where (∇*∇α)ᵢ is the component-wise rough Laplacian.

For form_degree=2 the operator is obtained from the scalar Laplace–Beltrami via the Hodge star:

Δ( f dx∧dy ) = (Δ₀( f / √|g| )) √|g| dx∧dy .

Parameters:

• metric (Metric) – Must have metric.dim == 2.

• form_degree ({0, 1, 2}) – Degree of the input differential form.

Returns:

• dict with keys

• * ``’principal’`` (sympy.Expr — leading symbol gⁱʲ ξᵢ ξⱼ.)

• * ``’subprincipal’`` (sympy.Expr — first-order transport term) – (non-zero for k=0; zero for k=1,2 at symbol level).

• * ``’full’`` (sympy.Expr — σ₂ + i σ₁ (microlocal convention).)

• * ``’weitzenbock’`` (sympy.Expr or None — the Weitzenböck curvature) – correction K (present for k=1,2; None for k=0).

• * ``’action’`` (callable — action(form) applies the operator) – symbolically and returns the result. - k=0: form is a scalar SymPy expression;

  returns a scalar Δf.

  • k=1: form is a 2‑tuple (α₀, α₁); returns a 2‑tuple ((Δα)₀, (Δα)₁).

  • k=2: form is a scalar SymPy expression f; returns a scalar representing the coefficient of Δ( f dx∧dy ) = (Δ( f/√g )) √g dx∧dy.

Raises:

NotImplementedError – If called on a 1D metric, or if form_degree is not 0, 1, or 2.

Examples

Flat metric — Laplacian of a 2‑form:

>>> from sympy import symbols, Matrix, simplify
>>> x, y = symbols('x y', real=True)
>>> m = Metric(Matrix([[1, 0], [0, 1]]), (x, y))
>>> op = de_rham_laplacian(m, form_degree=2)
>>> simplify(op['action'](x**2))   # Δ( x² dx∧dy ) = 2 dx∧dy
2

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

Compute the geodesic distance between two points on a 2D Riemannian manifold.

Two complementary numerical methods are provided:

• Shooting (default): iteratively refines an initial tangent vector at p by shooting a geodesic and comparing its endpoint to q. Converges quadratically near the solution and is accurate for well-separated points without conjugate loci between them.

• Optimisation: minimises the energy functional E(v) = ½ ‖v‖_g² + penalty · ‖exp_p(v) − q‖² over the initial tangent vector v using BFGS (scipy.optimize.minimize). More robust when the shooting iteration diverges, but generally less precise.

Parameters:

• metric (Metric) – Must have metric.dim == 2.

• p (tuple of two floats) – Start and end points (x, y) in the manifold.

• q (tuple of two floats) – Start and end points (x, y) in the manifold.

• method ({‘shooting’, ‘optimize’}, default 'shooting') – Numerical method to use.

• max_iter (int, default 50) – Maximum number of shooting iterations (ignored for 'optimize').

• tol (float, default 1e-6) – Convergence tolerance on the endpoint error ‖exp_p(v) − q‖ for the shooting method.

Returns:

Approximate geodesic distance d(p, q).

Return type:

float

Raises:

• NotImplementedError – If called on a 1D metric.

• ValueError – If method is not 'shooting' or 'optimize'.

Notes

Both methods rely on exponential_map() internally, so the accuracy of the geodesic integrator (controlled by n_steps inside exponential_map) directly affects the result. For very curved spaces or points near conjugate loci, neither method is guaranteed to converge.

Examples

Flat metric — Euclidean distance:

>>> import numpy as np
>>> x, y = symbols('x y', real=True)
>>> m = Metric(Matrix([[1, 0], [0, 1]]), (x, y))
>>> np.isclose(distance(m, (0, 0), (3, 4), method='shooting'), 5.0, rtol=1e-3)
True
>>> np.isclose(distance(m, (0, 0), (3, 4), method='optimize'), 5.0, rtol=5e-2)
True

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

Evaluate the Riemannian exponential map exp_p(t·v) on a 2D manifold.

The exponential map sends the tangent vector v at base point p along the unique unit-speed geodesic γ with γ(0) = p and γ̇(0) = v, returning the point γ(t). Geodesic completeness is assumed; the result may be inaccurate near conjugate loci.

Parameters:

• metric (Metric) – Must have metric.dim == 2.

• p (tuple of two floats) – Base point (x₀, y₀) in the manifold.

• v (tuple of two floats) – Initial tangent vector (vₓ, vᵧ) at p. The geodesic is parametrised so that ‖γ̇(0)‖_g = ‖v‖_g, i.e. the speed equals the norm of v.

• t (float, default 1.0) – Parameter value at which to evaluate the geodesic.

• method (str, default 'rk45') – Geodesic integration method passed to geodesic_solver().

Returns:

End point (x(t), y(t)) of the geodesic.

Return type:

tuple of two floats

Raises:

NotImplementedError – If called on a 1D metric.

Examples

Flat metric — exp_p(tv) = p + tv:

>>> x, y = symbols('x y', real=True)
>>> m = Metric(Matrix([[1, 0], [0, 1]]), (x, y))
>>> import numpy as np
>>> end = exponential_map(m, (0, 0), (3, 4), t=1.0)
>>> np.allclose(end, (3, 4), atol=1e-4)
True

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

Integrate geodesic flow using the Hamiltonian formulation with symplectic schemes.

Reformulates the geodesic equation as a Hamiltonian system with the kinetic-energy Hamiltonian

H = ½ gⁱʲ(q) pᵢ pⱼ,

where pᵢ = gᵢⱼ ẋʲ are the canonical momenta conjugate to the coordinates qⁱ. Hamilton’s equations

q̇ⁱ = ∂H/∂pᵢ = gⁱʲ pⱼ, ṗᵢ = −∂H/∂qⁱ

are then integrated by the companion symplectic.hamiltonian_flow module. Initial velocities are converted to momenta automatically.

Parameters:

• metric (Metric) – A 1D or 2D Riemannian metric.

• p0 (float (1D) or tuple of two floats (2D)) – Initial position in configuration space.

• v0 (float (1D) or tuple of two floats (2D)) – Initial velocity, converted internally to canonical momentum via pᵢ = gᵢⱼ(p0) v0ʲ.

• tspan (tuple of two floats) – Integration interval (t_start, t_end).

• method (str, default 'verlet') –

  Symplectic integration scheme forwarded to symplectic.hamiltonian_flow. The mapping from user-facing names to internal names is:

  • 'verlet' → Störmer–Verlet (2nd order, time-reversible).

  • 'stormer' → same as 'verlet' (alias).

  • 'symplectic' → same as 'verlet' (backward-compatible alias).

  • 'symplectic_euler' → symplectic Euler (1st order).

  • 'rk45' → Dormand–Prince (not symplectic, but high accuracy).

• n_steps (int, default 1000) – Number of integration steps.

Returns:

Phase-space trajectory. Keys by dimension:

• 1D: 't', 'x', 'v', 'p', 'energy'.

• 2D: 't', 'x', 'y', 'vx', 'vy', 'px', 'py', 'energy'.

The 'energy' array contains H evaluated at each time step; for a perfect symplectic integrator it is exactly conserved.

Return type:

dict

Raises:

• ValueError – If method is not in the allowed list.

• NotImplementedError – If metric.dim is neither 1 nor 2.

Notes

Symplectic integrators conserve a modified (shadow) Hamiltonian exactly, so energy will exhibit bounded oscillations rather than secular drift. The Störmer–Verlet scheme is second-order accurate and time-reversible, making it the default choice for long-time integrations.

Examples

1D energy conservation over 10 time units:

>>> x = symbols('x', real=True, positive=True)
>>> m = Metric(x**2, (x,))
>>> res = geodesic_hamiltonian_flow(m, 2.0, 2.5, (0, 10), method='verlet', n_steps=2000)
>>> import numpy as np
>>> np.std(res['energy']) / res['energy'][0] < 0.01   # < 1% drift
True

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

Integrate the geodesic equation on a 1D or 2D Riemannian manifold.

The geodesic equation in local coordinates is

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

supplemented by the initial conditions x(0) = p0, ẋ(0) = v0. This function dispatches automatically to _geodesic_1d or _geodesic_2d based on metric.dim.

Parameters:

• metric (Metric) – The Riemannian metric defining the manifold.

• p0 (float (1D) or tuple of two floats (2D)) – Initial position. For 2D pass (x₀, y₀).

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

• tspan (tuple of two floats) – Integration interval (t_start, t_end).

• method (str, default 'rk4') –

  Numerical integration scheme. Available methods by dimension:

  • 1D: 'rk4' (RK23 via SciPy), 'adaptive' (RK45), 'symplectic' (leapfrog/Störmer–Verlet in phase space).

  • 2D: 'rk45' (RK45), 'rk4' (RK23), 'symplectic' or 'verlet' (Störmer–Verlet via Hamiltonian flow).

• n_steps (int, default 1000) – Number of time steps (evaluation points for fixed-step methods; used as t_eval grid for adaptive methods).

• reparametrize (bool, default False) – 2D only. If True, also compute the accumulated arc-length parameter s(t) = ∫₀ᵗ ‖ẋ(τ)‖_g dτ and include it in the output dict under the key 'arc_length'.

Returns:

Trajectory arrays. Keys depend on dimension and method:

• 1D: 't', 'x', 'v'. Symplectic method additionally includes 'p' (canonical momentum).

• 2D: 't', 'x', 'y', 'vx', 'vy'. Hamiltonian methods additionally include 'px', 'py', 'energy'. If reparametrize=True, also includes 'arc_length'.

Return type:

dict

Raises:

ValueError – If method is not recognised for the given dimension.

Examples

1D flat metric — straight-line geodesic:

>>> import numpy as np
>>> x = symbols('x', real=True)
>>> m = Metric(1, (x,))
>>> traj = geodesic_solver(m, 0.0, 2.0, (0, 3.0))
>>> np.isclose(traj['x'][-1], 6.0, rtol=1e-4)
True

2D Euclidean metric with arc-length reparametrisation:

>>> x, y = symbols('x y', real=True)
>>> m = Metric(Matrix([[1, 0], [0, 1]]), (x, y))
>>> traj = geodesic_solver(m, (0, 0), (1, 1), (0, 5), reparametrize=True)
>>> 'arc_length' in traj
True

riemannian.hodge_decomposition(metric, omega_components, domain, resolution=50, form_degree=1)

Numerically decompose a differential k‑form into its Hodge components on a 2D rectangular domain.

Supported form degrees

form_degree=0

Decomposes a scalar function f (0‑form) as

f = Δu + h₀,

where Δu is the co‑exact part (the Laplacian of a scalar potential u) and h₀ is the constant harmonic part (the weighted mean of f). The exact part is zero (no (−1)-form). The potential u solves the Neumann problem Δu = f − h₀, with the gauge fixed by pinning one interior grid point to zero. The decomposition is then

coexact = Δu harmonic = h₀

The harmonic part corresponds to the de Rham cohomology class [f] ∈ H⁰_dR(M), which is simply the constant component.

form_degree=1 (default)

Decomposes a 1‑form α = αₓ dx + αᵧ dy as

α = dφ + ⋆dψ + h,

where dφ is exact, ⋆dψ is co‑exact, and h is harmonic (Δh = 0).

The scalar potentials φ and ψ satisfy the Poisson equations

Δ₀ φ = δα (Dirichlet boundary condition) Δ₀ ψ = -δ(⋆α) (Neumann boundary condition)

with the gauge fixed by pinning an interior grid point to zero for the Neumann problem. The decomposition is then

α_exact = ∇φ = dφ α_coexact = ⋆dψ α_harmonic = α − α_exact − α_coexact

The harmonic part corresponds to the de Rham cohomology class [α] ∈ H¹_dR(M).

form_degree=2

Decomposes a 2‑form ω = f dx∧dy into exact and harmonic parts:

ω = d(⋆dφ) + h,

(there is no co‑exact component because no 3‑forms exist in 2D). The algorithm solves a single Dirichlet problem for the scalar potential φ defined by Δ₀ φ = -⋆ω, where ⋆ω = f / √|g| is the Hodge dual (a 0‑form). The exact part is then

ω_exact = d(⋆dφ)

and the harmonic part is ω_harmonic = ω − ω_exact.

The resulting harmonic 2‑form belongs to the second de Rham cohomology class [ω] ∈ H²_dR(M).

param metric:

Must be 2D.

type metric:

Metric

param omega_components:

• form_degree=0: a scalar f (0‑form).

• form_degree=1: a 2‑tuple (αₓ, αᵧ).

• form_degree=2: a scalar f representing ω = f dx∧dy.

type omega_components:

tuple of two (callable or sympy.Expr), or scalar

param domain:

((x_min, x_max), (y_min, y_max))

type domain:

tuple

param resolution:

Number of grid points per axis.

type resolution:

int, default 50

param form_degree:

Degree of the input differential form.

type form_degree:

{0, 1, 2}, default 1

returns:

All cases contain a key 'grid' holding the RiemannianGrid instance used for the computation.

form_degree=0:

• 'potential_u' – array : potential u (Neumann, pinned).

• 'coexact' – array : Δu = f − h₀.

• 'harmonic' – array : constant harmonic part h₀.

form_degree=1:

• 'potential_phi' – array : φ potential (Dirichlet).

• 'potential_psi' – array : ψ potential (Neumann).

• 'alpha_exact' – tuple of two arrays : (dφ)ₓ, (dφ)ᵧ.

• 'alpha_coexact' – tuple of two arrays : (⋆dψ)ₓ, (⋆dψ)ᵧ.

• 'alpha_harmonic' – tuple of two arrays : hₓ, hᵧ.

form_degree=2:

• 'potential_phi' – array : φ such that Δ₀ φ = ⋆ω.

• 'omega_exact' – array : coefficient of the exact part.

• 'omega_harmonic' – array : coefficient of the harmonic part.

rtype:

dict

raises NotImplementedError:

If metric.dim != 2 or form_degree is not 0, 1, or 2.

Examples

0‑form — function with zero mean on the unit square:

>>> from sympy import symbols, Matrix
>>> import numpy as np
>>> x, y = symbols('x y', real=True)
>>> m = Metric(Matrix([[1, 0], [0, 1]]), (x, y))
>>> f = x**2 - y**2
>>> dec0 = hodge_decomposition(m, f, ((0, 1), (0, 1)), resolution=40,
...                            form_degree=0)
>>> # harmonic part ≈ 0, coexact part ≈ f

1‑form — rotation form on the flat torus (purely harmonic):

>>> dec1 = hodge_decomposition(m, (-y, x), ((0, 1), (0, 1)), resolution=40)
>>> # harmonic part carries ≈ 100% of L² energy

2‑form — constant vorticity ω = dx∧dy on the flat plane:

>>> dec2 = hodge_decomposition(m, 1, ((0, 1), (0, 1)),
...                            resolution=40, form_degree=2)
>>> # omega_exact ≈ 0, omega_harmonic ≈ 1 (contractible domain: b₂ = 0)

riemannian.hodge_decomposition_0form(grid, omega_components, _eval)

Decompose a scalar function (0‑form) into co‑exact and harmonic parts.

For a scalar field f = Δu + h₀, where:

• Δu is the co‑exact part (the Laplacian of a scalar potential u),

• h₀ is the constant harmonic part (the weighted mean of f).

The potential u solves the Neumann problem Δu = f - h₀, with the gauge fixed by pinning one interior grid point to zero.

Parameters:

• grid (RiemannianGrid) – The grid on which the decomposition is performed.

• omega_components (scalar (callable, sympy.Expr, or ndarray)) – The input 0‑form f.

• _eval (callable) – Helper that evaluates symbolic or callable inputs on the grid.

Returns:

• ‘potential_u’ : ndarray — scalar potential u.

• ’coexact’ : ndarray — Δu = f - h₀.

• ’harmonic’ : ndarray — constant harmonic part h₀.

• ’grid’ : RiemannianGrid — the grid used.

Return type:

dict

riemannian.hodge_decomposition_1form(grid, omega_components, _eval, _codf, _star1, dx, dy)

Decompose a 1‑form α = αₓ dx + αᵧ dy into exact, co‑exact, and harmonic parts.

The decomposition follows the Hodge–Helmholtz theorem on a compact manifold with boundary:

α = dφ + ⋆dψ + h,

where:

• dφ is exact,

• ⋆dψ is co‑exact,

• h is harmonic (Δh = 0).

The potentials φ and ψ are scalar functions satisfying the Poisson equations

Δ₀ φ = δα (exact potential) Δ₀ ψ = -δ(⋆α) (co‑exact potential)

with boundary conditions:

• φ satisfies homogeneous Dirichlet BC (φ = 0 on ∂M) to ensure uniqueness,

• ψ satisfies homogeneous Neumann BC (∂ₙψ = 0 on ∂M) for the same reason.

The gauge for ψ is fixed by pinning one interior grid point to zero.

Parameters:

• grid (RiemannianGrid) – The grid on which the decomposition is performed.

• omega_components (tuple of two (callable, sympy.Expr, or ndarray)) – The input 1‑form components (αₓ, αᵧ).

• _eval (callable) – Helper that evaluates symbolic or callable inputs on the grid.

• _codf (callable) – Numerical codifferential operator (δ) for 1‑forms.

• _star1 (callable) – Numerical Hodge star operator for 1‑forms.

• dx (float) – Grid spacings (used for finite‑differences).

• dy (float) – Grid spacings (used for finite‑differences).

Returns:

• ‘potential_phi’ : ndarray — scalar potential φ.

• ’potential_psi’ : ndarray — scalar potential ψ.

• ’alpha_exact’ : tuple of two ndarrays — (dφ)ₓ, (dφ)ᵧ.

• ’alpha_coexact’ : tuple of two ndarrays — (⋆dψ)ₓ, (⋆dψ)ᵧ.

• ’alpha_harmonic’ : tuple of two ndarrays — hₓ, hᵧ.

• ’grid’ : RiemannianGrid — the grid used.

Return type:

dict

riemannian.hodge_decomposition_2form(grid, omega_components, _eval, _codf, _star1, _star2, dx, dy)

Decompose a 2‑form ω = f dx∧dy into exact and harmonic parts.

On a 2‑dimensional oriented manifold with boundary, the Hodge decomposition of a 2‑form ω is

ω = d(⋆dφ) + h,

where:

• d(⋆dφ) is exact (the exterior derivative of a 1‑form),

• h is harmonic (Δh = 0).

The scalar potential φ satisfies the Poisson equation

Δ₀ φ = -⋆ω

with homogeneous Dirichlet boundary condition (φ = 0 on ∂M). This choice ensures a unique solution. The exact part is then computed as

ω_exact = d(⋆dφ),

and the harmonic part as ω_harmonic = ω − ω_exact.

Parameters:

• grid (RiemannianGrid) – The grid on which the decomposition is performed.

• omega_components (scalar (callable, sympy.Expr, or ndarray)) – The coefficient f of the 2‑form ω = f dx∧dy.

• _eval (callable) – Helper that evaluates symbolic or callable inputs on the grid.

• _codf (callable) – Numerical codifferential operator for 1‑forms (unused here but kept for signature consistency).

• _star1 (callable) – Numerical Hodge star operator for 1‑forms (used to obtain ⋆dφ).

• _star2 (callable) – Numerical Hodge star operator for 2‑forms (used to compute ⋆ω).

• dx (float) – Grid spacings (used for finite‑differences).

• dy (float) – Grid spacings (used for finite‑differences).

Returns:

• ‘potential_phi’ : ndarray — scalar potential φ.

• ’omega_exact’ : ndarray — coefficient of the exact part.

• ’omega_harmonic’ : ndarray — coefficient of the harmonic part.

• ’grid’ : RiemannianGrid — the grid used.

Return type:

dict

riemannian.hodge_star(metric, form_degree)

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

The Hodge star is the metric-induced isomorphism between k-forms and (n−k)-forms. In an oriented 2D Riemannian manifold with metric g and volume form dV = √|det g| dx∧dy, it acts as:

• 0-forms (functions): ⋆f = f √|det g|.

• 1-forms: if α = αₓ dx + αᵧ dy then ⋆α = (g⁻¹ antisymmetric contraction) · √|det g|, concretely

  ⋆α = (g⁰⁰ αᵧ − g⁰¹ αₓ) √|g| dx + (−g⁰¹ αᵧ + g¹¹ αₓ) √|g| dy.

• 2-forms: ⋆(f dx∧dy) = f / √|det g|.

Parameters:

• metric (Metric) – Must have metric.dim == 2.

• form_degree ({0, 1, 2}) – Degree k of the input differential form.

Returns:

• For degrees 0 and 2: a function of a single SymPy expression (the coefficient of the form) returning the coefficient of the Hodge dual.

• For degree 1: a function (alpha_x, alpha_y) → (beta_x, beta_y) returning the two components of the Hodge dual 1-form.

Return type:

callable

Raises:

• NotImplementedError – If called on a 1D metric.

• ValueError – If form_degree is not 0, 1, or 2.

Examples

Scaled Euclidean metric g = diag(4, 9):

>>> from sympy import symbols, Matrix, simplify
>>> x, y = symbols('x y', real=True)
>>> m = Metric(Matrix([[4, 0], [0, 9]]), (x, y))
>>> star0 = hodge_star(m, 0)
>>> simplify(star0(1))   # ⋆1 = √(4·9) = 6
6
>>> star2 = hodge_star(m, 2)
>>> simplify(star2(12))  # ⋆(12 dx∧dy) = 12/6 = 2
2

riemannian.induced_metric(R, du, dv)

Compute the metric tensor induced by a 3D embedding via finite differences.

Given a surface parameterised as R(u, v) ∈ ℝ³, the induced (pullback) metric is defined by

g_{ij} = ⟨∂_i R, ∂_j R⟩,

i.e. g11 = |∂R/∂u|², g12 = ⟨∂R/∂u, ∂R/∂v⟩, g22 = |∂R/∂v|².

Partial derivatives are approximated with second-order central differences in the interior and first-order one-sided differences at the boundary rows and columns.

Parameters:

• R (ndarray, shape (nu, nv, 3)) – Embedding array. R[i, j] is the position vector in ℝ³ at the grid point (u_i, v_j).

• du (float) – Uniform grid spacing in the u-direction.

• dv (float) – Uniform grid spacing in the v-direction.

Returns:

• g11 (ndarray, shape (nu, nv)) – ⟨∂R/∂u, ∂R/∂u⟩ at each grid point.

• g12 (ndarray, shape (nu, nv)) – ⟨∂R/∂u, ∂R/∂v⟩ at each grid point.

• g22 (ndarray, shape (nu, nv)) – ⟨∂R/∂v, ∂R/∂v⟩ at each grid point.

• dRdu (ndarray, shape (nu, nv, 3)) – Finite-difference approximation of ∂R/∂u.

• dRdv (ndarray, shape (nu, nv, 3)) – Finite-difference approximation of ∂R/∂v.

Notes

The returned dRdu and dRdv are reused by metric_deficit and corrugation_step to avoid redundant computation.

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

Solve the Jacobi equation (geodesic deviation equation) along a geodesic (2D only).

The Jacobi equation describes how a one-parameter family of geodesics spreads apart. For a reference geodesic γ(t) with tangent vector v = γ̇, the Jacobi field J satisfies

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

where R is the Riemann curvature tensor and D/dt denotes the covariant derivative along γ. In local coordinates (x⁰, x¹) this system is equivalent to the first‑order ODEs:

dJⁱ/dt = DJⁱ - Γⁱⱼₖ Jʲ vᵏ d(DJⁱ)/dt = -Rⁱⱼₖₗ Jʲ vᵏ vˡ - Γⁱⱼₖ DJʲ vᵏ

The geodesic trajectory (x(t), y(t), vx(t), vy(t)) is interpolated cubically from the input dict. The ODE is integrated with high‑accuracy BDF (backward differentiation formula) from scipy.integrate.solve_ivp.

Parameters:

• metric (Metric) – Must have metric.dim == 2.

• geodesic (dict) – Output dict from geodesic_solver() with keys 't', 'x', 'y', 'vx', 'vy'. The time array must span at least tspan.

• initial_variation (dict) –

  Initial conditions for the Jacobi field:

  • 'J0' : tuple of two floats — initial value (Jₓ(0), Jᵧ(0)).

  • 'DJ0' : tuple of two floats — initial covariant derivative (DJₓ/dt(0), DJᵧ/dt(0)).

• tspan (tuple of two floats) – Integration interval for the Jacobi ODE. Must be a sub‑interval of the time range covered by geodesic['t'].

• n_steps (int, default 1000) – Number of output time steps (equally spaced in the interval tspan).

Returns:

Keys:

• 't' : ndarray — time values.

• 'J_x', 'J_y' : ndarray — Jacobi field components.

• 'DJ_x', 'DJ_y' : ndarray — covariant derivative of the Jacobi field.

Return type:

dict

Raises:

NotImplementedError – If called on a 1D metric.

Notes

Conjugate points along the geodesic occur where the Jacobi field vanishes. Large growth indicates negative sectional curvature in the geodesic direction, while oscillatory behaviour indicates positive curvature.

The ODE is solved with high precision (rtol=1e-9, atol=1e-12) and a stiff solver (BDF) to handle possible rapid changes in curvature.

In 2D the Riemann tensor has only one independent component. All curvature terms are evaluated via the Gauss equation R(J,v)v = K(g(v,v)J - g(J,v)v), where K is computed numerically along the geodesic via the Brioschi formula on a local grid. This replaces the symbolic riemann_tensor() call and eliminates its lambdify overhead.

Examples

>>> import numpy as np
>>> from sympy import symbols, Matrix, sin
>>> theta, phi = symbols('theta phi', real=True)
>>> m = Metric(Matrix([[1, 0], [0, sin(theta)**2]]), (theta, phi))
>>> geod = geodesic_solver(m, (np.pi/2, 0), (0, 1), (0, 2), n_steps=200)
>>> jac = jacobi_equation_solver(m, geod, {'J0': (0, 0), 'DJ0': (0.1, 0)}, (0, 2))
>>> 'J_x' in jac and 'J_y' in jac
True

riemannian.laplace_beltrami(metric)

Return the symbol dict of the Laplace–Beltrami operator for a metric.

Convenience wrapper around Metric.laplace_beltrami_symbol(). The returned dict contains the principal symbol (the inverse-metric quadratic form in the cotangent variables), the subprincipal symbol (the transport term), and their complex combination as the full microlocal symbol.

Parameters:

metric (Metric) – A 1D or 2D Riemannian metric.

Returns:

With keys 'principal', 'subprincipal', and 'full'. See Metric.laplace_beltrami_symbol() for full documentation.

Return type:

dict

Examples

>>> x = symbols('x', real=True, positive=True)
>>> m = Metric(x**2, (x,))
>>> lb = laplace_beltrami(m)
>>> lb['principal']   # g¹¹ ξ² = ξ²/x²
xi**2/x**2

riemannian.metric_deficit(R, g11_t, g12_t, g22_t, du, dv)

Compute the pointwise metric deficit and its Frobenius norm.

The deficit is the difference between the target metric (the abstract Riemannian tensor we wish to realise) and the metric actually induced by the current embedding R:

Δg = g_target − g_induced(R).

The scalar summary is the root-mean-square Frobenius norm

‖Δg‖_F = sqrt( mean( Δg11² + 2·Δg12² + Δg22² ) ),

where the factor of 2 on the off-diagonal term accounts for symmetry. This quantity monotonically decreasing toward zero is the convergence criterion used by add_corrugations.

Parameters:

• R (ndarray, shape (nu, nv, 3)) – Current embedding.

• g11_t (ndarray, shape (nu, nv)) – Target metric component g_{uu}.

• g12_t (ndarray, shape (nu, nv)) – Target metric component g_{uv}.

• g22_t (ndarray, shape (nu, nv)) – Target metric component g_{vv}.

• du (float) – Uniform grid spacing in the u-direction.

• dv (float) – Uniform grid spacing in the v-direction.

Returns:

• dg11 (ndarray, shape (nu, nv)) – Pointwise deficit in the g_{uu} component.

• dg12 (ndarray, shape (nu, nv)) – Pointwise deficit in the g_{uv} component.

• dg22 (ndarray, shape (nu, nv)) – Pointwise deficit in the g_{vv} component.

• frob (float) – RMS Frobenius norm of the deficit matrix over the grid.

riemannian.parallel_transport(metric, curve, initial_vector, tspan=None, method='RK45')

Transport a vector along a curve using parallel transport.

The parallel transport equation Dv/dt = 0 is solved as a linear ODE:

dv^i/dt = - Γ^i_{jk} v^j ẋ^k,

where ẋ^k are the components of the curve’s velocity.

When the curve dict was produced by geodesic_solver() it already contains exact velocity arrays ('v' in 1D, 'vx' / 'vy' in 2D). These are used directly, avoiding the numerical differentiation noise that np.gradient would introduce. For curves supplied without velocity arrays the function falls back to np.gradient automatically.

Parameters:

• metric (Metric) – The Riemannian metric.

• curve (dict) – A trajectory dict from geodesic_solver containing at least: - ‘t’ : array of time values - for 1D: ‘x’ array - for 2D: ‘x’, ‘y’ arrays

• initial_vector (float (1D) or tuple of two floats (2D)) – Components of the vector to transport at t = curve[‘t’][0].

• tspan (tuple (t0, t1) or None) – Time interval over which to transport. If None, the whole curve is used.

• method (str, default 'RK45') – Integration method passed to solve_ivp.

Returns:

• ‘t’ : array of times (same as curve[‘t’] within tspan)

• For 1D: ‘v’ : transported vector components

• For 2D: ‘vx’, ‘vy’ : transported vector components

Return type:

dict

riemannian.plot_corrugation_pipeline(result, title='', colormap='plasma', dark=True)

Visualise the full Nash–Kuiper corrugation pipeline in a single figure.

The figure has two rows:

• Top row: Up to 5 evenly spaced 3D surface snapshots from the pipeline (short map → intermediate iterations → final embedding), each labelled with its iteration index and Frobenius deficit.

• Bottom row: A line plot of the Frobenius metric deficit versus iteration number, showing convergence toward zero.

Parameters:

• result (dict) – The dictionary returned by add_corrugations. Must contain keys 'snapshots' (list of ndarray) and 'deficits' (list of float).

• title (str, optional) – Overall figure title (suptitle). Default is no title.

• colormap (str, optional) – Any Matplotlib colormap name applied to z-height colouring. Default is 'plasma'.

• dark (bool, optional) – If True (default), dark background with white text. If False, white background with black text.

Returns:

fig – The figure object. Its size adapts to the number of snapshots shown: width = 4 × n_snaps inches, height = 9 inches.

Return type:

matplotlib.figure.Figure

Notes

At most 5 snapshots are shown regardless of how many iterations were run. The indices are chosen by np.linspace so the first (short map) and last (final) snapshots are always included.

riemannian.plot_embedding(R, title='', colormap='plasma', dark=True, backend='widget')

Render a single 3D embedding as a shaded surface coloured by z-height.

The z-values are linearly normalised to [0, 1] before being mapped through the chosen colormap, so the full color range is always used regardless of the absolute scale of the embedding.

Parameters:

• R (ndarray, shape (nu, nv, 3)) – Embedding to display. Axes 0 and 1 are the u- and v-parameter directions; axis 2 holds the (x, y, z) coordinates.

• title (str, optional) – Figure title displayed above the axes. Default is no title.

• colormap (str, optional) – Any Matplotlib colormap name. Default is 'plasma'.

• dark (bool, optional) – If True (default), use a near-black background and white labels, suitable for publication or presentation slides. If False, use a white background with black labels.

• backend (str, optional) – Matplotlib backend to use for this function only. Options: ‘widget’ (interactive), ‘inline’ (static), or None (use global backend).

Returns:

• fig (matplotlib.figure.Figure) – The figure object.

• ax (matplotlib.axes.Axes3D) – The 3D axes.

riemannian.sturm_liouville_reduce(metric, potential_expr=None)

Reduce the Laplace–Beltrami eigenvalue problem to Sturm–Liouville form (1D only).

The Laplace–Beltrami eigenvalue problem on a 1D Riemannian manifold,

−Δ_g u + V u = λ u,

is equivalent to the classical Sturm–Liouville problem

−(p u′)′ + q u = λ w u,

with weight function w = √g, coefficient p = √g · g¹¹, and potential coefficient q = V √g. This standard form enables use of classical spectral theory, finite-element methods, and Sturm–Liouville solvers.

Parameters:

• metric (Metric) – Must have metric.dim == 1.

• potential_expr (sympy.Expr or None, default None) – Optional potential V(x). Pass None (default) for the pure Laplace–Beltrami operator (V = 0).

Returns:

• dict with keys

• * ``’p’`` (sympy.Expr — coefficient p(x) = √g · g¹¹.)

• * ``’q’`` (sympy.Expr — coefficient q(x) = V(x) √g (zero if no potential).)

• * ``’w’`` (sympy.Expr — weight function w(x) = √g.)

• * ``’p_func’`` (callable — numerical p(x_val).)

• * ``’q_func’`` (callable — numerical q(x_val).)

• * ``’w_func’`` (callable — numerical w(x_val).)

Raises:

NotImplementedError – If called on a 2D metric.

Examples

Cone metric g = x², no potential:

>>> x = symbols('x', real=True, positive=True)
>>> m = Metric(x**2, (x,))
>>> sl = sturm_liouville_reduce(m)
>>> simplify(sl['p'])   # p = √(x²) · (1/x²) = 1/x
1/x
>>> simplify(sl['w'])   # w = √(x²) = x
x
>>> sl['q']             # q = 0 (no potential)
0

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

Numerically verify the Gauss–Bonnet theorem over a rectangular domain (2D only).

The Gauss–Bonnet theorem states that for a compact oriented surface M without boundary,

∫∫_M K dA = 2π χ(M),

where K is the Gaussian curvature, dA = √|det g| dx dy is the Riemannian area element, and χ(M) is the Euler characteristic. For the unit sphere χ = 2, giving ∫K dA = 4π; for a torus χ = 0, giving ∫K dA = 0.

This function numerically evaluates the left-hand side over a rectangular domain using scipy.integrate.dblquad and compares it to the expected value 2π (assuming χ = 1 for the supplied region, e.g. a topological disk).

Parameters:

• metric (Metric) – Must have metric.dim == 2.

• domain (tuple) – ((x_min, x_max), (y_min, y_max)) — rectangular integration region.

• resolution (int, default 100) – Unused by the current implementation (kept for API consistency; the adaptive dblquad routine controls its own resolution internally).

Returns:

• dict with keys

• * ``’integral’`` (float — numerically computed ∫∫_M K dA.)

• * ``’expected’`` (float — reference value 2π.)

• * ``’integration_error’`` (float — absolute error estimate from dblquad.)

• * ``’relative_error’`` (float — |integral − 2π| / 2π.)

Raises:

NotImplementedError – If called on a 1D metric.

Notes

Accuracy degrades near coordinate singularities (e.g. θ = 0 or π on the sphere) where the curvature expression or the metric determinant may become ill-conditioned numerically. Restrict the domain away from such points.

Examples

>>> x, y = symbols('x y', real=True)
>>> m = Metric(Matrix([[1, 0], [0, 1]]), (x, y))    # flat: K = 0
>>> result = verify_gauss_bonnet(m, ((0, 1), (0, 1)))
>>> abs(result['integral']) < 1e-10   # ∫K dA = 0 for flat metric
True

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

Visualise a curvature-related quantity of a 1D or 2D Riemannian manifold.

Dispatches to _visualize_curvature_1d or _visualize_curvature_2d based on metric.dim and renders a Matplotlib figure.

Parameters:

• metric (Metric) – The Riemannian metric to visualise.

• x_range (tuple of two floats, optional) – Coordinate range along the first axis for the plot. For 1D, inferred from geodesic extents if initial_conditions is also passed; defaults to (-5, 5) if neither is provided. For 2D, required.

• y_range (tuple of two floats, optional) – Coordinate range along the second axis. Required for 2D metrics. Unused for 1D.

• resolution (int, default 100) – Number of sample points along each axis when building the background curve (1D) or colour mesh (2D).

• quantity (str, default 'gauss') –

  Which scalar field to display:

  • 2D: 'gauss' — Gaussian curvature K(x, y); 'ricci_scalar' — Ricci scalar R(x, y).

  • 1D: 'metric' or 'gauss' — metric component g₁₁(x); 'christoffel' — Christoffel symbol Γ¹₁₁(x).

• cmap (str, default 'RdBu_r') – Matplotlib colormap name for the 2D colour mesh. Ignored for 1D.

• **kwargs –

  Additional keyword arguments forwarded to the 1D renderer only:

  • initial_conditions : list of (x₀, v₀) — if provided, geodesics are overlaid on the lower subplot.

  • tspan : tuple — integration interval for geodesics (default (0, 10)).

  • colorby : str — colouring scheme for geodesics (default 'speed').

  • n_steps : int — number of integration steps (default 500).

Returns:

Displays a Matplotlib figure via plt.show().

Return type:

None

Raises:

ValueError – If x_range or y_range are missing for a 2D metric, or if quantity is not recognised for the given dimension.

Examples

2D — Gaussian curvature of the Poincaré half-plane:

>>> x, y = symbols('x y', real=True)
>>> m = Metric(Matrix([[1/y**2, 0], [0, 1/y**2]]), (x, y))
>>> visualize_curvature(m, x_range=(-2, 2), y_range=(0.5, 3), quantity='gauss')

1D — metric and geodesics for g = x²:

>>> x = symbols('x', real=True, positive=True)
>>> m = Metric(x**2, (x,))
>>> visualize_curvature(m, x_range=(0.5, 5), quantity='metric',
...                     initial_conditions=[(1.0, 0.5)], tspan=(0, 4))

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

Visualise geodesic trajectories on a 1D or 2D Riemannian manifold.

Dispatches to _visualize_geodesics_1d or _visualize_geodesics_2d based on metric.dim. Each geodesic is integrated from the given initial conditions and plotted using Matplotlib.

Parameters:

• metric (Metric) – The Riemannian metric defining the manifold.

• initial_conditions (list) –

  • 1D: list of (x₀, v₀) pairs (floats).

  • 2D: list of ((x₀, y₀), (vₓ₀, vᵧ₀)) pairs of tuples.

• tspan (tuple of two floats) – Integration interval (t_start, t_end) for all geodesics.

• x_range (tuple of two floats, optional) – Horizontal plotting range. Inferred from trajectory extents if None.

• y_range (tuple of two floats, optional) – Vertical plotting range (2D only). Inferred from trajectories if None.

• colorby ({‘speed’, ‘time’, ‘curvature’}, default 'speed') – 1D only. Quantity used to colour the scatter plot of each geodesic point. 'speed' uses |ẋ|, 'time' uses t, and 'curvature' uses |Γ¹₁₁(x(t))|.

• plot_curvature (bool, default True) – 2D only. If True, render the Gaussian curvature as a colour-mapped background heatmap behind the geodesic curves.

• n_steps (int, default 500) – Number of integration time steps per geodesic.

Returns:

Displays a Matplotlib figure via plt.show().

Return type:

None

Examples

1D — visualise two geodesics on the cone metric g = x²:

>>> x = symbols('x', real=True, positive=True)
>>> m = Metric(x**2, (x,))
>>> visualize_geodesics(m, [(1.0, 0.5), (2.0, -1.0)], (0, 5))

2D — geodesics on the Poincaré half-plane:

>>> x, y = symbols('x y', real=True)
>>> m = Metric(Matrix([[1/y**2, 0], [0, 1/y**2]]), (x, y))
>>> visualize_geodesics(m, [((0, 1), (0, 1)), ((0, 1), (1, 0))], (0, 2))

riemannian.visualize_hodge_decomposition(decomp, domain=None, resolution=50, cmap='RdBu_r', quiver_stride=3, form_degree=None, use_3d=False)

Visualise the Hodge decomposition of a differential form on a 2D manifold.

This function works for decompositions of 0‑, 1‑, and 2‑forms as returned by hodge_decomposition(). It automatically detects the form degree from the keys present in decomp, or you can specify it explicitly with the form_degree parameter.

For a 0‑form (scalar function) f = Δu + h₀, the visualisation shows:

• Top row: original scalar field, co‑exact part Δu, harmonic part h₀.

• Bottom row: potential u, residual error, and an energy bar chart comparing the two contributions.

For a 1‑form α = dφ + ⋆dψ + h, the visualisation shows:

• Top row: original vector field and its three components (exact, co‑exact, harmonic) as quiver plots overlaid on a heatmap of the magnitude.

• Bottom row: scalar potentials φ and ψ, residual error, and an energy bar chart.

For a 2‑form ω = d(⋆dφ) + h (no co‑exact component), the display consists of:

• Top row: original scalar field, exact part, harmonic part.

• Bottom row: scalar potential φ, residual error, and an energy bar chart comparing the exact and harmonic contributions.

The energy distribution uses the metric‑weighted L² norm:

• For 1‑forms: ∫ ‖α‖_g² dV

• For 0‑forms and 2‑forms: ∫ f² dV

This requires that the decomposition dictionary contains a ‘grid’ key with a RiemannianGrid instance that provides the metric weights. If the grid is missing, the unweighted sum of squared values is used.

Parameters:

• decomp (dict) –

  Output of hodge_decomposition(). Must contain keys appropriate for the form degree:

  • 0‑form: ‘potential_u’, ‘coexact’, ‘harmonic’.

  • 1‑form: ‘alpha_exact’, ‘alpha_coexact’, ‘alpha_harmonic’, ‘potential_phi’, ‘potential_psi’.

  • 2‑form: ‘omega_exact’, ‘omega_harmonic’, ‘potential_phi’.

• domain (tuple, optional) – ((x_min, x_max), (y_min, y_max)) – coordinate bounds. If the decomposition does not contain a ‘grid’ key, this parameter is required to build the coordinate mesh. If a grid is present, domain is ignored.

• resolution (int, default 50) – Number of grid points per axis. Used only when a grid is not available in decomp.

• cmap (str, default 'RdBu_r') – Matplotlib colormap for scalar fields (and magnitude for 1‑forms).

• quiver_stride (int, default 3) – Sub‑sampling stride for quiver arrows (only used for 1‑forms).

• form_degree ({0, 1, 2}, optional) – Explicitly specify the form degree. If None, the degree is inferred from the dictionary keys.

Returns:

Displays a Matplotlib figure.

Return type:

None

Examples

0‑form on the unit square:

>>> from sympy import symbols, Matrix
>>> x, y = symbols('x y', real=True)
>>> m = Metric(Matrix([[1, 0], [0, 1]]), (x, y))
>>> f = x**2 - y**2
>>> dec0 = hodge_decomposition(m, f, ((0, 1), (0, 1)), form_degree=0)
>>> visualize_hodge_decomposition(dec0)

1‑form on the flat torus:

>>> dec1 = hodge_decomposition(m, (-y, x), ((0, 2*np.pi), (0, 2*np.pi)))
>>> visualize_hodge_decomposition(dec1)

2‑form on the unit square (constant vorticity):

>>> dec2 = hodge_decomposition(m, 1, ((0, 1), (0, 1)), form_degree=2)
>>> visualize_hodge_decomposition(dec2)
use_3d : bool, default False
    If True and a grid is available, render on a 3D embedding of the
    surface.

Returns:

Displays a Matplotlib figure.

Return type:

None