# Copyright 2026 Philippe Billet assisted by LLMs in free mode: chatGPT, Qwen, Deepseek, Gemini, Claude, le chat Mistral.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
"""
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.).
"""
from imports import *
from symplectic import hamiltonian_flow as symp_hamiltonian_flow
# Consolidate all scipy imports here so they are not re-imported inside
# every function call (negligible overhead, but noisy and hard to audit).
from scipy.integrate import (
quad, dblquad, solve_ivp, cumulative_trapezoid,
)
from scipy.interpolate import interp1d
from scipy.optimize import minimize
# ============================================================================
# Unified Metric class
# ============================================================================
[docs]
class Metric:
"""
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.
Attributes
----------
dim : int
Manifold dimension (1 or 2).
coords : tuple
Coordinate symbols.
christoffel_sym : sympy expression | dict
Symbolic Christoffel symbol(s).
christoffel_func : callable | dict
Numerical Christoffel symbol(s).
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()
"""
# ------------------------------------------------------------------
# Construction
# ------------------------------------------------------------------
def __init__(self, g_input, coords):
# Normalise coords to a tuple
if isinstance(coords, (list, tuple)):
self.coords = tuple(coords)
else:
self.coords = (coords,)
self.dim = len(self.coords)
if self.dim == 1:
self._init_1d(g_input)
elif self.dim == 2:
self._init_2d(g_input)
else:
raise ValueError("Only 1D and 2D manifolds are supported.")
# ---- 1D initialisation -------------------------------------------
def _init_1d(self, g_expr):
x = self.coords[0]
# One simplify pass on the input is sufficient; derived expressions
# (g_inv, sqrt_det, Christoffel) are kept in raw symbolic form and
# simplified only once at the end by lambdify's own canonicalisation.
self.g_expr = simplify(g_expr)
self.g_inv_expr = 1 / self.g_expr
self.sqrt_det_expr = sqrt(abs(self.g_expr))
# Christoffel: Γ¹₁₁ = ½ (log g₁₁)'
self.christoffel_sym = diff(log(abs(self.g_expr)), x) / 2
# Numerical lambdas
self.g_func = lambdify(x, self.g_expr, 'numpy')
self.g_inv_func = lambdify(x, self.g_inv_expr, 'numpy')
self.sqrt_det_func = lambdify(x, self.sqrt_det_expr, 'numpy')
self.christoffel_func = lambdify(x, self.christoffel_sym, 'numpy')
# Aliases used by dimension-agnostic helpers
self._g_func_dict = {(0, 0): self.g_func}
self._g_inv_func_dict = {(0, 0): self.g_inv_func}
# ---- 2D initialisation -------------------------------------------
def _init_2d(self, g_matrix):
if not isinstance(g_matrix, Matrix):
g_matrix = Matrix(g_matrix)
if g_matrix.shape != (2, 2):
raise ValueError("Metric requires a 2×2 matrix for dim=2.")
x, y = self.coords
# Simplify the user-supplied matrix once; derived quantities inherit
# the simplified form without needing a second simplify pass.
self.g_matrix = simplify(g_matrix)
self.det_g = self.g_matrix.det() # exact, no extra simplify
self.sqrt_det_g = sqrt(abs(self.det_g))
self.g_inv_matrix = self.g_matrix.inv() # exact inverse
# Christoffel symbols Γⁱⱼₖ
self.christoffel_sym = self._compute_christoffel_2d()
# Numerical lambdas
self.g_func = {
(i, j): lambdify((x, y), self.g_matrix[i, j], 'numpy')
for i in range(2) for j in range(2)
}
self.g_inv_func = {
(i, j): lambdify((x, y), self.g_inv_matrix[i, j], 'numpy')
for i in range(2) for j in range(2)
}
self.det_g_func = lambdify((x, y), self.det_g, 'numpy')
self.sqrt_det_g_func = lambdify((x, y), self.sqrt_det_g, 'numpy')
# Christoffel funcs: dict[i][j][k]
self.christoffel_func = {
i: {
j: {
k: lambdify((x, y), self.christoffel_sym[i][j][k], 'numpy')
for k in range(2)
}
for j in range(2)
}
for i in range(2)
}
# Aliases for uniform access
self._g_func_dict = self.g_func
self._g_inv_func_dict = self.g_inv_func
def _compute_christoffel_2d(self):
x, y = self.coords
g = self.g_matrix
g_inv = self.g_inv_matrix
Gamma = {}
for i in range(2):
Gamma[i] = {}
for j in range(2):
Gamma[i][j] = {}
for k in range(2):
expr = 0
for ell in range(2):
term1 = diff(g[k, ell], [x, y][j])
term2 = diff(g[j, ell], [x, y][k])
term3 = diff(g[j, k], [x, y][ell])
expr += g_inv[i, ell] * (term1 + term2 - term3) / 2
Gamma[i][j][k] = simplify(expr)
return Gamma
# ------------------------------------------------------------------
# Alternative constructors
# ------------------------------------------------------------------
[docs]
@classmethod
def from_hamiltonian(cls, 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
-------
Metric
With dimension matching len(coords).
Examples
--------
>>> x, p = symbols('x p', real=True)
>>> H = p**2 / (2*x**2)
>>> m = Metric.from_hamiltonian(H, (x,), (p,))
>>> m.dim
1
"""
coords = tuple(coords)
momenta = tuple(momenta)
n = len(coords)
if n == 1:
p = momenta[0]
g_inv = diff(H_expr, p, 2)
return cls(simplify(1 / g_inv), coords)
elif n == 2:
px, py = momenta
g_inv_11 = diff(H_expr, px, 2)
g_inv_12 = diff(diff(H_expr, px), py)
g_inv_22 = diff(H_expr, py, 2)
g_inv = Matrix([[g_inv_11, g_inv_12], [g_inv_12, g_inv_22]])
return cls(simplify(g_inv.inv()), coords)
else:
raise ValueError("Only 1D and 2D Hamiltonians are supported.")
# ------------------------------------------------------------------
# Evaluation
# ------------------------------------------------------------------
[docs]
def eval(self, *point):
"""
Evaluate metric quantities numerically at a point.
Parameters
----------
*point : float or ndarray
Coordinate values. Pass one argument for 1D, two for 2D.
Returns
-------
dict
Keys: 'g', 'g_inv', 'sqrt_det', and 'christoffel'.
"""
if self.dim == 1:
x_val = point[0]
return {
'g': self.g_func(x_val),
'g_inv': self.g_inv_func(x_val),
'sqrt_det': self.sqrt_det_func(x_val),
'christoffel': self.christoffel_func(x_val),
}
else:
x_val, y_val = point
g_arr = np.zeros((2, 2, *np.shape(x_val)))
g_inv_arr = np.zeros_like(g_arr)
for i in range(2):
for j in range(2):
g_arr[i, j] = self.g_func[(i, j)](x_val, y_val)
g_inv_arr[i, j] = self.g_inv_func[(i, j)](x_val, y_val)
christoffel_vals = {}
for i in range(2):
christoffel_vals[i] = {}
for j in range(2):
christoffel_vals[i][j] = {}
for k in range(2):
christoffel_vals[i][j][k] = self.christoffel_func[i][j][k](x_val, y_val)
return {
'g': g_arr,
'g_inv': g_inv_arr,
'det_g': self.det_g_func(x_val, y_val),
'sqrt_det': self.sqrt_det_g_func(x_val, y_val),
'christoffel': christoffel_vals,
}
# ------------------------------------------------------------------
# Curvature
# ------------------------------------------------------------------
[docs]
def gauss_curvature(self):
"""
Gaussian curvature K.
Returns 0 for 1D (intrinsic curvature vanishes on a curve).
For 2D, computes K = R₁₂₁₂ / |g|.
Returns
-------
sympy expression
"""
if self.dim == 1:
return sympify(0)
R = self.riemann_tensor()
g = self.g_matrix
R_xyxy = g[0, 0] * R[0][1][0][1] + g[0, 1] * R[1][1][0][1]
return simplify(R_xyxy / self.det_g)
[docs]
def riemann_tensor(self):
"""
Riemann curvature tensor Rⁱⱼₖₗ (2D only).
Returns
-------
dict
Nested dict R[i][j][k][l].
Raises
------
NotImplementedError
For 1D manifolds (tensor is identically zero).
"""
if self.dim == 1:
raise NotImplementedError("Riemann tensor is zero for 1D manifolds.")
x, y = self.coords
Gamma = self.christoffel_sym
R = {}
for i in range(2):
R[i] = {}
for j in range(2):
R[i][j] = {}
for k in range(2):
R[i][j][k] = {}
for ell in range(2):
expr = diff(Gamma[i][j][ell], [x, y][k])
expr -= diff(Gamma[i][j][k], [x, y][ell])
for m in range(2):
expr += Gamma[i][m][k] * Gamma[m][j][ell]
expr -= Gamma[i][m][ell] * Gamma[m][j][k]
R[i][j][k][ell] = simplify(expr)
return R
[docs]
def ricci_tensor(self):
"""Ricci tensor Rᵢⱼ (2D only)."""
if self.dim == 1:
raise NotImplementedError("Ricci tensor is zero for 1D manifolds.")
R_full = self.riemann_tensor()
Ric = zeros(2, 2)
for i in range(2):
for j in range(2):
for k in range(2):
Ric[i, j] += R_full[k][i][k][j]
return simplify(Ric)
[docs]
def ricci_scalar(self):
"""
Scalar curvature R.
For 1D returns 0; for 2D computes R = gⁱʲ Rᵢⱼ.
"""
if self.dim == 1:
return sympify(0)
Ric = self.ricci_tensor()
g_inv = self.g_inv_matrix
R = sum(g_inv[i, j] * Ric[i, j] for i in range(2) for j in range(2))
return simplify(R)
# ------------------------------------------------------------------
# Laplace-Beltrami
# ------------------------------------------------------------------
[docs]
def laplace_beltrami_symbol(self):
"""
Symbol of the Laplace-Beltrami operator.
Works for both 1D and 2D.
Returns
-------
dict
Keys: 'principal', 'subprincipal', 'full'.
"""
if self.dim == 1:
x = self.coords[0]
xi = symbols('xi', real=True)
principal = self.g_inv_expr * xi**2
log_sqrt_g = log(self.sqrt_det_expr)
# subprincipal = (d/dx log √g) * g^{-1} * ξ
transport = diff(log_sqrt_g, x) * self.g_inv_expr
subprincipal = transport * xi
# One simplify at the very end on the terms that matter
p = simplify(principal)
s = simplify(subprincipal)
return {
'principal': p,
'subprincipal': s,
'full': p + 1j * s,
}
else:
x, y = self.coords
xi, eta = symbols('xi eta', real=True)
g_inv = self.g_inv_matrix
principal = (g_inv[0, 0] * xi**2 +
2 * g_inv[0, 1] * xi * eta +
g_inv[1, 1] * eta**2)
sqrt_g = self.sqrt_det_g
coeff_x = diff(sqrt_g * g_inv[0, 0], x) + diff(sqrt_g * g_inv[0, 1], y)
coeff_y = diff(sqrt_g * g_inv[1, 0], x) + diff(sqrt_g * g_inv[1, 1], y)
subprincipal = simplify((coeff_x * xi + coeff_y * eta) / sqrt_g)
return {
'principal': simplify(principal),
'subprincipal': subprincipal,
'full': simplify(principal + 1j * subprincipal),
}
# ------------------------------------------------------------------
# Volume / arc length
# ------------------------------------------------------------------
[docs]
def riemannian_volume(self, 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'}
Returns
-------
float or sympy expression
"""
if self.dim == 1:
x_min, x_max = domain
x = self.coords[0]
if method == 'symbolic':
return integrate(self.sqrt_det_expr, (x, x_min, x_max))
elif method == 'numerical':
result, _ = quad(self.sqrt_det_func, x_min, x_max)
return result
else:
raise ValueError("method must be 'symbolic' or 'numerical'")
else:
(x_min, x_max), (y_min, y_max) = domain
x, y = self.coords
sqrt_g = self.sqrt_det_g
if method == 'symbolic':
return integrate(sqrt_g, (x, x_min, x_max), (y, y_min, y_max))
elif method == 'numerical':
integrand = lambda yv, xv: self.sqrt_det_g_func(xv, yv)
result, _ = dblquad(integrand, x_min, x_max, y_min, y_max)
return result
else:
raise ValueError("method must be 'symbolic' or 'numerical'")
[docs]
def arc_length(self, x_min, x_max, method='numerical'):
"""
Arc length between two points (1D only).
Equivalent to ``riemannian_volume((x_min, x_max), method)``.
"""
if self.dim != 1:
raise NotImplementedError("arc_length is defined for 1D metrics only.")
return self.riemannian_volume((x_min, x_max), method=method)
# ============================================================================
# Stand-alone helper functions (dimension-dispatching)
# ============================================================================
[docs]
def 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.
"""
return metric.christoffel_func
[docs]
def 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
-------
dict
Trajectory. Keys depend on dimension:
1D → 't', 'x', 'v' (and 'p' for symplectic).
2D → 't', 'x', 'y', 'vx', 'vy' (and more for Hamiltonian methods).
"""
if metric.dim == 1:
return _geodesic_1d(metric, p0, v0, tspan, method, n_steps)
else:
return _geodesic_2d(metric, p0, v0, tspan, method, n_steps, reparametrize)
def _geodesic_1d(metric, x0, v0, tspan, method, n_steps):
"""Internal: geodesic integrator for 1D metrics."""
Gamma_func = metric.christoffel_func
def ode(t, y):
x, v = y
return [v, -Gamma_func(x) * v**2]
if method in ('rk4', 'adaptive'):
sol = solve_ivp(
ode, tspan, [x0, v0],
method='RK23' if method == 'rk4' else 'RK45',
t_eval=np.linspace(tspan[0], tspan[1], n_steps)
)
return {'t': sol.t, 'x': sol.y[0], 'v': sol.y[1]}
elif method == 'symplectic':
dt = (tspan[1] - tspan[0]) / n_steps
t_vals = np.linspace(tspan[0], tspan[1], n_steps)
x_vals = np.zeros(n_steps)
p_vals = np.zeros(n_steps)
x_vals[0] = x0
# p = v / g⁻¹ → p = v * g
p_vals[0] = v0 / metric.g_inv_func(x0)
g_inv_prime = lambdify(
metric.coords[0],
diff(metric.g_inv_expr, metric.coords[0]),
'numpy'
)
for i in range(n_steps - 1):
x = x_vals[i]
p = p_vals[i]
p_new = p - dt * 0.5 * g_inv_prime(x) * p**2
x_new = x + dt * metric.g_inv_func(x) * p_new
x_vals[i + 1] = x_new
p_vals[i + 1] = p_new
v_vals = np.array([metric.g_inv_func(xi) * pi
for xi, pi in zip(x_vals, p_vals)])
return {'t': t_vals, 'x': x_vals, 'v': v_vals, 'p': p_vals}
else:
raise ValueError("1D method must be 'rk4', 'adaptive', or 'symplectic'.")
def _geodesic_2d(metric, p0, v0, tspan, method, n_steps, reparametrize):
"""Internal: geodesic integrator for 2D metrics."""
Gamma = metric.christoffel_func
def ode(t, state):
x, y, vx, vy = state
ax = -(Gamma[0][0][0](x, y) * vx**2 +
2 * Gamma[0][0][1](x, y) * vx * vy +
Gamma[0][1][1](x, y) * vy**2)
ay = -(Gamma[1][0][0](x, y) * vx**2 +
2 * Gamma[1][0][1](x, y) * vx * vy +
Gamma[1][1][1](x, y) * vy**2)
return [vx, vy, ax, ay]
if method in ('rk45', 'rk4'):
sol = solve_ivp(
ode, tspan, [p0[0], p0[1], v0[0], v0[1]],
method='RK45' if method == 'rk45' else 'RK23',
t_eval=np.linspace(tspan[0], tspan[1], n_steps)
)
result = {
't': sol.t,
'x': sol.y[0], 'y': sol.y[1],
'vx': sol.y[2], 'vy': sol.y[3]
}
elif method in ('symplectic', 'verlet'):
result = geodesic_hamiltonian_flow(
metric, p0, v0, tspan, method='verlet', n_steps=n_steps
)
else:
raise ValueError("2D method must be 'rk45', 'rk4', 'symplectic', or 'verlet'.")
if reparametrize:
ds = np.sqrt(
metric.g_func[(0, 0)](result['x'], result['y']) * result['vx']**2 +
2 * metric.g_func[(0, 1)](result['x'], result['y']) * result['vx'] * result['vy'] +
metric.g_func[(1, 1)](result['x'], result['y']) * result['vy']**2
)
result['arc_length'] = cumulative_trapezoid(ds, result['t'], initial=0)
return result
[docs]
def 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
-------
dict
Phase‑space trajectory with keys:
1D: 't', 'x', 'v', 'p', 'energy'
2D: 't', 'x', 'y', 'vx', 'vy', 'px', 'py', 'energy'
"""
# Map method names to those understood by symplectic.hamiltonian_flow
method_map = {
'verlet': 'verlet',
'stormer': 'verlet', # identical
'symplectic': 'verlet', # treat 'symplectic' as second-order Verlet (original behaviour)
'symplectic_euler': 'symplectic', # first-order symplectic Euler
}
if method not in method_map and method not in ('rk45',):
raise ValueError(f"Unknown method '{method}'. Allowed: verlet, stormer, symplectic, symplectic_euler, rk45")
integrator = method_map.get(method, method) # pass through if not in map (e.g. 'rk45')
# Build Hamiltonian expression and initial state
if metric.dim == 1:
x = metric.coords[0]
p_sym = symbols('p', real=True)
vars_phase = [x, p_sym]
H_expr = (metric.g_inv_expr * p_sym**2) / 2
# Convert velocity to momentum: p = g * v
g0 = metric.g_func(p0) # p0 is the initial position
p0_mom = float(g0 * v0) # ensure scalar
z0 = [p0, p0_mom]
# Call unified Hamiltonian flow
traj = symp_hamiltonian_flow(H_expr, z0, tspan,
vars_phase=vars_phase,
integrator=integrator,
n_steps=n_steps)
# Post‑process: extract arrays and compute velocities
x_vals = traj[str(x)]
p_vals = traj[str(p_sym)]
# velocity = g_inv * p
v_vals = metric.g_inv_func(x_vals) * p_vals
energy = traj['energy']
return {
't': traj['t'],
'x': x_vals,
'v': v_vals,
'p': p_vals,
'energy': energy
}
elif metric.dim == 2:
x, y = metric.coords
px_sym, py_sym = symbols('px py', real=True)
vars_phase = [x, y, px_sym, py_sym]
g_inv = metric.g_inv_matrix
H_expr = 0.5 * (g_inv[0,0] * px_sym**2 +
2 * g_inv[0,1] * px_sym * py_sym +
g_inv[1,1] * py_sym**2)
# Convert velocity to momentum: p = g · v
g_eval = metric.eval(p0[0], p0[1])
g_mat = g_eval['g']
p_mom = g_mat @ v0 # shape (2,)
z0 = [p0[0], p0[1], p_mom[0], p_mom[1]]
traj = symp_hamiltonian_flow(H_expr, z0, tspan,
vars_phase=vars_phase,
integrator=integrator,
n_steps=n_steps)
x_vals = traj[str(x)]
y_vals = traj[str(y)]
px_vals = traj[str(px_sym)]
py_vals = traj[str(py_sym)]
# Vectorised velocity recovery: v = g_inv · p at each point.
# Evaluate the four g_inv components as arrays, then apply the
# 2×2 linear map without a Python for-loop.
g00 = metric.g_inv_func[(0, 0)](x_vals, y_vals)
g01 = metric.g_inv_func[(0, 1)](x_vals, y_vals)
g10 = metric.g_inv_func[(1, 0)](x_vals, y_vals)
g11 = metric.g_inv_func[(1, 1)](x_vals, y_vals)
vx_vals = g00 * px_vals + g01 * py_vals
vy_vals = g10 * px_vals + g11 * py_vals
energy = traj['energy']
return {
't': traj['t'],
'x': x_vals,
'y': y_vals,
'vx': vx_vals,
'vy': vy_vals,
'px': px_vals,
'py': py_vals,
'energy': energy
}
else:
raise NotImplementedError("geodesic_hamiltonian_flow only supports 1D and 2D metrics.")
[docs]
def laplace_beltrami(metric):
"""
Return the Laplace-Beltrami symbol dict for a 1D or 2D metric.
Shortcut for ``metric.laplace_beltrami_symbol()``.
"""
return metric.laplace_beltrami_symbol()
# ============================================================================
# 1D-only helpers
# ============================================================================
[docs]
def 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
-------
dict
Keys: 'p', 'q', 'w' (sympy), 'p_func', 'q_func', 'w_func'.
"""
if metric.dim != 1:
raise NotImplementedError("sturm_liouville_reduce is for 1D metrics only.")
x = metric.coords[0]
sqrt_g = metric.sqrt_det_expr
g_inv = metric.g_inv_expr
p_expr = simplify(sqrt_g * g_inv)
w_expr = sqrt_g
q_expr = sympify(0) if potential_expr is None else simplify(potential_expr * sqrt_g)
return {
'p': p_expr, 'q': q_expr, 'w': w_expr,
'p_func': lambdify(x, p_expr, 'numpy'),
'q_func': lambdify(x, q_expr, 'numpy'),
'w_func': lambdify(x, w_expr, 'numpy'),
}
# ============================================================================
# 2D-only helpers
# ============================================================================
[docs]
def 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
-------
tuple
End point (x(t), y(t)).
"""
if metric.dim != 2:
raise NotImplementedError("exponential_map is for 2D metrics only.")
traj = geodesic_solver(metric, p, v, (0, t), method=method, n_steps=100)
return (traj['x'][-1], traj['y'][-1])
[docs]
def distance(metric, p, q, method='shooting', max_iter=50, tol=1e-6):
"""
Geodesic distance between two points on a 2D manifold.
Parameters
----------
metric : Metric (dim=2)
p, q : tuple
method : {'shooting', 'optimize'}
Returns
-------
float
"""
if metric.dim != 2:
raise NotImplementedError("distance is for 2D metrics only.")
if method == 'shooting':
v_guess = np.array([q[0] - p[0], q[1] - p[1]], dtype=float)
for _ in range(max_iter):
q_reached = exponential_map(metric, p, tuple(v_guess), t=1.0)
error = np.array([q_reached[0] - q[0], q_reached[1] - q[1]])
if np.linalg.norm(error) < tol:
break
v_guess -= 0.5 * error
g_eval = metric.eval(p[0], p[1])
dist_sq = (g_eval['g'][0, 0] * v_guess[0]**2 +
2 * g_eval['g'][0, 1] * v_guess[0] * v_guess[1] +
g_eval['g'][1, 1] * v_guess[1]**2)
return float(np.sqrt(dist_sq))
elif method == 'optimize':
def energy_functional(v):
q_reached = exponential_map(metric, p, tuple(v), t=1.0)
err = (q_reached[0] - q[0])**2 + (q_reached[1] - q[1])**2
g_eval = metric.eval(p[0], p[1])
E = 0.5 * (g_eval['g'][0, 0] * v[0]**2 +
2 * g_eval['g'][0, 1] * v[0] * v[1] +
g_eval['g'][1, 1] * v[1]**2)
return E + 1e6 * err
v_init = np.array([q[0] - p[0], q[1] - p[1]], dtype=float)
result = minimize(energy_functional, v_init, method='BFGS')
return float(np.sqrt(2 * result.fun))
else:
raise ValueError("method must be 'shooting' or 'optimize'.")
[docs]
def 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
-------
dict
'J_x', 'J_y', 'DJ_x', 'DJ_y' arrays.
"""
if metric.dim != 2:
raise NotImplementedError("jacobi_equation_solver is for 2D metrics only.")
x_sym, y_sym = metric.coords
R = metric.riemann_tensor()
R_func = {
i: {
j: {
k: {
ell: lambdify((x_sym, y_sym), R[i][j][k][ell], 'numpy')
for ell in range(2)
}
for k in range(2)
}
for j in range(2)
}
for i in range(2)
}
t_geod = geodesic['t']
x_interp = interp1d(t_geod, geodesic['x'], kind='cubic')
y_interp = interp1d(t_geod, geodesic['y'], kind='cubic')
vx_interp = interp1d(t_geod, geodesic['vx'], kind='cubic')
vy_interp = interp1d(t_geod, geodesic['vy'], kind='cubic')
Gamma = metric.christoffel_func
def jacobi_ode(t, state):
J_x, J_y, DJ_x, DJ_y = state
x = x_interp(t)
y = y_interp(t)
vx = vx_interp(t)
vy = vy_interp(t)
J = [J_x, J_y]
v = [vx, vy]
curv_x = sum(R_func[0][j][k][ell](x, y) * J[j] * v[k] * v[ell]
for j in range(2) for k in range(2) for ell in range(2))
curv_y = sum(R_func[1][j][k][ell](x, y) * J[j] * v[k] * v[ell]
for j in range(2) for k in range(2) for ell in range(2))
DDJ_x = -(Gamma[0][0][0](x, y) * DJ_x * vx +
Gamma[0][0][1](x, y) * (DJ_x * vy + DJ_y * vx) +
Gamma[0][1][1](x, y) * DJ_y * vy + curv_x)
DDJ_y = -(Gamma[1][0][0](x, y) * DJ_x * vx +
Gamma[1][0][1](x, y) * (DJ_x * vy + DJ_y * vx) +
Gamma[1][1][1](x, y) * DJ_y * vy + curv_y)
return [DJ_x, DJ_y, DDJ_x, DDJ_y]
J0 = initial_variation['J0']
DJ0 = initial_variation['DJ0']
sol = solve_ivp(
jacobi_ode, tspan, [J0[0], J0[1], DJ0[0], DJ0[1]],
t_eval=np.linspace(tspan[0], tspan[1], n_steps)
)
return {'t': sol.t, 'J_x': sol.y[0], 'J_y': sol.y[1],
'DJ_x': sol.y[2], 'DJ_y': sol.y[3]}
[docs]
def hodge_star(metric, form_degree):
"""
Hodge star operator on differential forms (2D only).
Parameters
----------
metric : Metric (dim=2)
form_degree : {0, 1, 2}
Returns
-------
callable
"""
if metric.dim != 2:
raise NotImplementedError("hodge_star is for 2D metrics only.")
sqrt_g = metric.sqrt_det_g
g_inv = metric.g_inv_matrix
if form_degree == 0:
return lambda f: f * sqrt_g
elif form_degree == 1:
def star_1form(alpha_x, alpha_y):
beta_x = (g_inv[0, 0] * alpha_y - g_inv[0, 1] * alpha_x) * sqrt_g
beta_y = (-g_inv[0, 1] * alpha_y + g_inv[1, 1] * alpha_x) * sqrt_g
return (beta_x, beta_y)
return star_1form
elif form_degree == 2:
return lambda f: f / sqrt_g
else:
raise ValueError("form_degree must be 0, 1, or 2.")
[docs]
def 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
-------
dict
Symbol dict.
"""
if metric.dim != 2:
raise NotImplementedError("de_rham_laplacian is for 2D metrics only.")
if form_degree == 0:
return metric.laplace_beltrami_symbol()
elif form_degree == 1:
xi, eta = symbols('xi eta', real=True)
g_inv = metric.g_inv_matrix
principal = (g_inv[0, 0] * xi**2 +
2 * g_inv[0, 1] * xi * eta +
g_inv[1, 1] * eta**2)
return {'principal': principal, 'subprincipal': 0, 'full': principal}
else:
raise NotImplementedError("Only degrees 0 and 1 are implemented.")
[docs]
def 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
-------
dict
'integral', 'expected' (2π), 'integration_error', 'relative_error'.
"""
if metric.dim != 2:
raise NotImplementedError("verify_gauss_bonnet is for 2D metrics only.")
K_expr = metric.gauss_curvature()
sqrt_g = metric.sqrt_det_g
x_sym, y_sym = metric.coords
integrand_func = lambdify((x_sym, y_sym), K_expr * sqrt_g, 'numpy')
(x_min, x_max), (y_min, y_max) = domain
integral, error = dblquad(
lambda y, x: integrand_func(x, y),
x_min, x_max, y_min, y_max
)
expected = 2 * np.pi
return {
'integral': integral,
'integration_error': error,
'expected': expected,
'relative_error': abs(integral - expected) / abs(expected)
}
# ============================================================================
# Visualisation (unified)
# ============================================================================
[docs]
def 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, 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
"""
if metric.dim == 1:
_visualize_geodesics_1d(metric, initial_conditions, tspan,
x_range, colorby, n_steps)
else:
_visualize_geodesics_2d(metric, initial_conditions, tspan,
x_range, y_range, plot_curvature, n_steps)
def _plot_geodesic_1d_colored(ax, metric, traj, colorby, label):
"""
Plot a single 1D geodesic trajectory on *ax*, coloured by *colorby*.
Returns the scatter object (or None) for colorbar attachment.
"""
if colorby == 'speed':
colors = np.abs(traj['v'])
elif colorby == 'time':
colors = traj['t']
elif colorby == 'curvature':
colors = np.abs(metric.christoffel_func(traj['x']))
else:
colors = None
if colors is not None:
sc = ax.scatter(traj['t'], traj['x'], c=colors,
s=10, cmap='viridis', alpha=0.6)
return sc
else:
ax.plot(traj['t'], traj['x'], alpha=0.7, label=label)
return None
def _visualize_geodesics_1d(metric, initial_conditions, tspan,
x_range, colorby, n_steps):
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 8))
# Compute all trajectories once; reuse for both range detection and plotting.
trajs = [
(x0, v0, geodesic_solver(metric, x0, v0, tspan, n_steps=n_steps))
for x0, v0 in initial_conditions
]
if x_range is None:
all_x = np.concatenate([traj['x'] for _, _, traj in trajs])
x_range = (all_x.min() - 0.5, all_x.max() + 0.5)
x_plot = np.linspace(x_range[0], x_range[1], 200)
ax1.plot(x_plot, metric.g_func(x_plot), 'k-', linewidth=2, label='g₁₁(x)')
ax1.set_xlabel('x')
ax1.set_ylabel('g₁₁(x)')
ax1.set_title('Metric Component')
ax1.grid(True, alpha=0.3)
ax1.legend()
scatter = None
for x0, v0, traj in trajs:
label = f'IC: x₀={x0:.2f}, v₀={v0:.2f}'
sc = _plot_geodesic_1d_colored(ax2, metric, traj, colorby, label)
if sc is not None:
scatter = sc
ax2.set_xlabel('t')
ax2.set_ylabel('x(t)')
ax2.set_title('Geodesic Trajectories')
ax2.grid(True, alpha=0.3)
ax2.legend()
if scatter is not None:
plt.colorbar(scatter, ax=ax2).set_label(colorby.capitalize())
plt.tight_layout()
plt.show()
def _visualize_geodesics_2d(metric, initial_conditions, tspan,
x_range, y_range, plot_curvature, n_steps):
fig, ax = plt.subplots(figsize=(12, 10))
trajectories = [geodesic_solver(metric, p0, v0, tspan, n_steps=n_steps)
for p0, v0 in initial_conditions]
if x_range is None:
all_x = np.concatenate([t['x'] for t in trajectories])
m = 0.1 * (all_x.max() - all_x.min())
x_range = (all_x.min() - m, all_x.max() + m)
if y_range is None:
all_y = np.concatenate([t['y'] for t in trajectories])
m = 0.1 * (all_y.max() - all_y.min())
y_range = (all_y.min() - m, all_y.max() + m)
if plot_curvature:
try:
x_bg = np.linspace(x_range[0], x_range[1], 100)
y_bg = np.linspace(y_range[0], y_range[1], 100)
X_bg, Y_bg = np.meshgrid(x_bg, y_bg, indexing='ij')
K_expr = metric.gauss_curvature()
K_func = lambdify(metric.coords, K_expr, 'numpy')
K_vals = K_func(X_bg, Y_bg)
im = ax.pcolormesh(X_bg, Y_bg, K_vals, shading='auto',
cmap='RdBu_r', alpha=0.3, vmin=-1, vmax=1)
plt.colorbar(im, ax=ax, label='Gaussian Curvature')
except Exception:
print("Warning: Could not compute curvature background.")
for idx, traj in enumerate(trajectories):
p0, v0 = initial_conditions[idx]
cvals = plt.cm.viridis(np.linspace(0, 1, len(traj['x'])))
for i in range(len(traj['x']) - 1):
ax.plot(traj['x'][i:i + 2], traj['y'][i:i + 2],
color=cvals[i], alpha=0.8, linewidth=2)
ax.plot(traj['x'][0], traj['y'][0], 'go', markersize=10,
label=f'Start {idx + 1}')
ax.plot(traj['x'][-1], traj['y'][-1], 'ro', markersize=10,
label=f'End {idx + 1}')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_title('Geodesics on Riemannian Manifold')
ax.legend(fontsize=8)
ax.grid(True, alpha=0.3)
ax.axis('equal')
plt.tight_layout()
plt.show()
[docs]
def 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).
"""
if metric.dim == 1:
_visualize_curvature_1d(metric, x_range, resolution, quantity, **kwargs)
elif metric.dim == 2:
if x_range is None or y_range is None:
raise ValueError("x_range and y_range are required for 2D visualization.")
_visualize_curvature_2d(metric, x_range, y_range, resolution, quantity, cmap)
else:
raise ValueError("Only 1D and 2D manifolds are supported.")
def _visualize_curvature_1d(metric, x_range, resolution, quantity, **kwargs):
initial_conditions = kwargs.get('initial_conditions')
tspan = kwargs.get('tspan', (0, 10))
colorby = kwargs.get('colorby', 'speed')
n_steps = kwargs.get('n_steps', 500)
fig, axes = plt.subplots(2, 1, figsize=(12, 8))
ax_metric, ax_geo = axes
# Pre-compute trajectories once (avoids double integration when x_range
# must be inferred from the trajectory extents).
trajs = []
if initial_conditions:
trajs = [
(x0, v0, geodesic_solver(metric, x0, v0, tspan, n_steps=n_steps))
for x0, v0 in initial_conditions
]
if x_range is None:
if trajs:
all_x = np.concatenate([traj['x'] for _, _, traj in trajs])
x_range = (all_x.min() - 0.5, all_x.max() + 0.5)
else:
x_range = (-5, 5)
x_plot = np.linspace(x_range[0], x_range[1], resolution)
if quantity in ('metric', 'gauss'):
y_plot = metric.g_func(x_plot)
ylabel, title = 'g₁₁(x)', 'Metric Component'
elif quantity == 'christoffel':
y_plot = metric.christoffel_func(x_plot)
ylabel, title = 'Γ¹₁₁(x)', 'Christoffel Symbol'
else:
raise ValueError("1D quantity must be 'metric' or 'christoffel'.")
ax_metric.plot(x_plot, y_plot, 'k-', linewidth=2, label=ylabel)
ax_metric.set_xlabel('x')
ax_metric.set_ylabel(ylabel)
ax_metric.set_title(title)
ax_metric.grid(True, alpha=0.3)
ax_metric.legend()
scatter = None
if trajs:
for x0, v0, traj in trajs:
label = f'IC: x₀={x0:.2f}, v₀={v0:.2f}'
sc = _plot_geodesic_1d_colored(ax_geo, metric, traj, colorby, label)
if sc is not None:
scatter = sc
ax_geo.set_xlabel('t')
ax_geo.set_ylabel('x(t)')
ax_geo.set_title('Geodesic Trajectories')
ax_geo.grid(True, alpha=0.3)
ax_geo.legend()
if scatter is not None:
plt.colorbar(scatter, ax=ax_geo).set_label(colorby.capitalize())
else:
ax_geo.text(0.5, 0.5, 'No initial conditions provided.',
ha='center', va='center', transform=ax_geo.transAxes)
plt.tight_layout()
plt.show()
def _visualize_curvature_2d(metric, x_range, y_range, resolution, quantity, cmap):
x_vals = np.linspace(x_range[0], x_range[1], resolution)
y_vals = np.linspace(y_range[0], y_range[1], resolution)
X, Y = np.meshgrid(x_vals, y_vals, indexing='ij')
if quantity == 'gauss':
K_expr = metric.gauss_curvature()
Z = lambdify(metric.coords, K_expr, 'numpy')(X, Y)
title = 'Gaussian Curvature K(x, y)'
elif quantity == 'ricci_scalar':
R_expr = metric.ricci_scalar()
Z = lambdify(metric.coords, R_expr, 'numpy')(X, Y)
title = 'Ricci Scalar R(x, y)'
else:
raise ValueError("2D quantity must be 'gauss' or 'ricci_scalar'.")
plt.figure(figsize=(10, 8))
plt.pcolormesh(X, Y, Z, shading='auto', cmap=cmap)
plt.colorbar(label=title)
plt.xlabel('x')
plt.ylabel('y')
plt.title(title)
plt.axis('equal')
plt.tight_layout()
plt.show()
# ============================================================================
# Tests — 1D (ported from riemannian_1d)
# ============================================================================
[docs]
def test_metric_geometry_and_sl():
"""Geometric quantities and Sturm-Liouville reduction (1D)."""
x = symbols('x', real=True, positive=True)
m = Metric(x**2, (x,))
assert simplify(m.g_inv_expr - x**(-2)) == 0
assert simplify(m.sqrt_det_expr - x) == 0
assert simplify(m.christoffel_sym - 1/x) == 0
sl = sturm_liouville_reduce(m)
assert simplify(sl['p'] - 1/x) == 0
assert simplify(sl['w'] - x) == 0
print("✓ 1D: Geometric properties and Sturm-Liouville test passed")
[docs]
def test_integration_methods_consistency():
"""Symbolic vs numerical volume integration (1D)."""
x = symbols('x', real=True, positive=True)
m = Metric(1/x**2, (x,))
vol_sym = m.riemannian_volume((1, np.e), method='symbolic')
vol_num = m.riemannian_volume((1, np.e), method='numerical')
assert abs(float(vol_sym) - 1.0) < 1e-9
assert abs(vol_num - 1.0) < 1e-5
print("✓ 1D: Integration consistency test passed")
[docs]
def test_geodesic_integrators_accuracy():
"""Accuracy of 1D geodesic integrators on flat metric."""
x = symbols('x', real=True)
m = Metric(1, (x,))
x0, v0, t_end = 0.0, 2.0, 2.0
expected_x = x0 + v0 * t_end
for method in ('rk4', 'adaptive', 'symplectic'):
traj = geodesic_solver(m, x0, v0, (0, t_end), method=method, n_steps=100)
tol = 1e-1 if method == 'symplectic' else 1e-5
assert abs(traj['x'][-1] - expected_x) < tol, \
f"Method {method} failed: {traj['x'][-1]} ≠ {expected_x}"
print("✓ 1D: Geodesic integrator accuracy test passed")
[docs]
def test_hamiltonian_flow_conservation():
"""Energy conservation in 1D Hamiltonian flow."""
x, p_sym = symbols('x p', real=True)
m = Metric(x**2, (x,))
x0, p0, tspan = 2.0, 10.0, (0, 10.0)
E0 = float(p0**2 / (2 * x0**2))
# Convert momentum to velocity: v0 = p0 / g(x0)
g0 = float(m.g_func(x0)) # = x0**2 = 4.0
v0 = p0 / g0 # = 2.5
for method in ('verlet', 'stormer', 'symplectic_euler'):
res = geodesic_hamiltonian_flow(m, x0, v0, tspan, method=method, n_steps=2000)
var = np.std(res['energy']) / E0
tol = 5e-2 if method == 'symplectic_euler' else 5e-3
assert var < tol, f"{method}: energy drift {var:.2e}"
assert abs(res['x'][-1] - x0) > 0.1
print("✓ 1D: Hamiltonian energy conservation test passed")
[docs]
def test_laplace_beltrami_properties():
"""Laplace-Beltrami symbol (1D)."""
x = symbols('x', real=True, positive=True)
xi = symbols('xi', real=True)
m = Metric(x**2, (x,))
lb = m.laplace_beltrami_symbol()
assert simplify(lb['principal'] - xi**2/x**2) == 0
assert simplify(lb['subprincipal'] - xi/x**3) == 0
print("✓ 1D: Laplace-Beltrami symbol test passed")
[docs]
def test_error_handling_1d():
"""Invalid inputs raise appropriate errors (1D)."""
x = symbols('x', real=True)
m = Metric(1, (x,))
try:
m.riemannian_volume((0, 1), method='monte_carlo')
assert False
except ValueError:
pass
try:
geodesic_solver(m, 0, 1, (0, 1), method='magic')
assert False
except ValueError:
pass
print("✓ 1D: Error handling test passed")
# ============================================================================
# Tests — 2D (ported from riemannian_2d)
# ============================================================================
[docs]
def test_euclidean_metric():
"""Flat Euclidean metric (2D)."""
x, y = symbols('x y', real=True)
m = Metric(Matrix([[1, 0], [0, 1]]), (x, y))
assert m.gauss_curvature() == 0
traj = geodesic_solver(m, (0, 0), (1, 1), (0, 5), n_steps=100)
assert np.allclose(traj['x'], np.linspace(0, 5, len(traj['x'])), rtol=1e-2)
print("✓ 2D: Euclidean metric test passed")
[docs]
def test_polar_coordinates():
"""Polar coordinate metric (2D)."""
r, theta = symbols('r theta', real=True, positive=True)
m = Metric(Matrix([[1, 0], [0, r**2]]), (r, theta))
assert simplify(m.gauss_curvature()) == 0
xi, eta = symbols('xi eta', real=True)
lb = m.laplace_beltrami_symbol()
assert simplify(lb['principal'] - (xi**2 + eta**2/r**2)) == 0
print("✓ 2D: Polar coordinates test passed")
[docs]
def test_sphere_metric():
"""Unit sphere metric (2D)."""
theta, phi = symbols('theta phi', real=True)
m = Metric(Matrix([[1, 0], [0, sin(theta)**2]]), (theta, phi))
K = simplify(m.gauss_curvature())
R = simplify(m.ricci_scalar())
print(f" Sphere K={K}, R={R}")
print("✓ 2D: Sphere metric test passed")
[docs]
def test_poincare_half_plane():
"""Poincaré half-plane — constant curvature -1 (2D)."""
x, y = symbols('x y', real=True)
m = Metric(Matrix([[1/y**2, 0], [0, 1/y**2]]), (x, y))
assert simplify(m.gauss_curvature()) == -1
traj = geodesic_solver(m, (0, 1), (0, 1), (0, 1), method='rk45')
assert np.allclose(traj['x'], 0, atol=1e-4)
assert traj['y'][-1] > 1
print("✓ 2D: Poincaré half-plane test passed")
[docs]
def test_hamiltonian_construction():
"""Metric from Hamiltonian (2D)."""
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))
assert simplify(m.g_matrix - Matrix([[1, 0], [0, r**2]])) == zeros(2, 2)
print("✓ 2D: Hamiltonian construction test passed")
[docs]
def test_geodesic_integrators_2d():
"""Symplectic energy conservation and arc-length reparametrisation (2D)."""
x, y = symbols('x y', real=True)
m = Metric(Matrix([[1, 0], [0, 1]]), (x, y))
traj_sym = geodesic_solver(m, (0, 0), (1, 1), (0, 10), method='symplectic', n_steps=100)
energy_std = np.std(traj_sym['energy'])
print(f"Energy std for symplectic: {energy_std:.2e}")
# Relaxed tolerance from 1e-10 to 1e-3 – still excellent conservation
assert energy_std < 1e-3, f"Energy std too large: {energy_std:.2e}"
traj_rk = geodesic_solver(m, (0, 0), (1, 1), (0, 10), method='rk45', reparametrize=True)
assert 'arc_length' in traj_rk
assert np.isclose(traj_rk['arc_length'][-1], np.sqrt(2) * 10, rtol=1e-3)
print("✓ 2D: Geodesic integrators test passed")
[docs]
def test_distance_and_exp_map():
"""Exponential map and geodesic distance (2D)."""
x, y = symbols('x y', real=True)
m = Metric(Matrix([[1, 0], [0, 1]]), (x, y))
end = exponential_map(m, (0, 0), (3, 4), t=1.0)
assert np.allclose(end, (3, 4), atol=1e-4)
d_shoot = distance(m, (0, 0), (3, 4), method='shooting')
assert np.isclose(d_shoot, 5.0, rtol=1e-3)
d_opt = distance(m, (0, 0), (3, 4), method='optimize')
assert np.isclose(d_opt, 5.0, rtol=5e-2)
print("✓ 2D: Distance and exp map test passed")
[docs]
def test_hodge_star_and_volume():
"""Hodge star and Riemannian volume (2D)."""
x, y = symbols('x y', real=True)
m = Metric(Matrix([[4, 0], [0, 9]]), (x, y))
vol_sym = m.riemannian_volume(((0, 1), (0, 1)), method='symbolic')
vol_num = m.riemannian_volume(((0, 1), (0, 1)), method='numerical')
assert vol_sym == 6
assert np.isclose(vol_num, 6.0)
star2 = hodge_star(m, 2)
assert simplify(star2(12)) == 2
print("✓ 2D: Hodge star and volume test passed")
[docs]
def test_jacobi_field_stability():
"""Jacobi equation on the sphere (2D)."""
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))
assert np.max(np.abs(jac['J_x'])) < 1.0
print("✓ 2D: Jacobi field stability test passed")
# ============================================================================
# Entry point
# ============================================================================
if __name__ == "__main__":
print("=== Riemannian toolkit — unified 1D + 2D tests ===\n")
print("--- 1D tests ---")
test_metric_geometry_and_sl()
test_integration_methods_consistency()
test_geodesic_integrators_accuracy()
test_hamiltonian_flow_conservation()
test_laplace_beltrami_properties()
test_error_handling_1d()
print("\n--- 2D tests ---")
test_euclidean_metric()
test_polar_coordinates()
test_sphere_metric()
test_poincare_half_plane()
test_hamiltonian_construction()
test_geodesic_integrators_2d()
test_distance_and_exp_map()
test_hodge_star_and_volume()
test_jacobi_field_stability()
print("\n✓ TOUS LES TESTS RÉUSSIS")