∑MathAdvanced

Smooth Manifolds & Tangent Spaces

Key Points

•
A smooth manifold is a space that looks like ordinary Euclidean space when you zoom in, glued together using charts that transition smoothly.
•
An atlas is a collection of charts whose overlap maps (transition maps) are smooth; this smoothness is the essence of the manifold’s differentiable structure.
•
The tangent space at a point is the set of all possible directions you can move from that point within the manifold, forming a vector space of the same dimension as the manifold.
•
Tangent vectors can be defined equivalently as velocities of curves, as derivations acting on smooth functions, or via Jacobians of charts.
•
In coordinates, tangent vectors transform by multiplying with the Jacobian of the transition map; this is how “components” change between charts.
•
For embedded manifolds like the sphere $S^{2}$ ⊂ $R^{3}$ , tangent spaces can also be seen as vectors orthogonal to the normal (e.g., v ⋅ $p = 0$ on $S^{2}$ ).
•
Computationally, you can implement charts, transition maps, and Jacobians in C++ to transform tangent vectors and verify smooth compatibility numerically.
•
Stereographic projection provides concrete charts for $S^{2}$ , with closed-form transition maps that make excellent test cases for code.

Prerequisites

→Multivariable Calculus (Derivatives and Jacobians) — Charts and transition maps use partial derivatives and the chain rule; Jacobians are essential for coordinate changes.
→Linear Algebra (Vector Spaces and Matrices) — Tangent spaces are vector spaces; Jacobians, projections, and least-squares require matrix operations.
→Topology (Basic Notions) — Manifold definitions rely on open sets, continuity, and homeomorphisms.
→Ordinary Differential Equations (Curves and Velocities) — Viewing tangent vectors as curve velocities needs understanding of differentiable curves.
→Numerical Methods (Finite Differences) — Computing Jacobians in code often uses finite-difference approximations and stability heuristics.

Detailed Explanation

Tap terms for definitions

01Overview

Hook: Imagine navigating a curved planet with flat maps. No single map covers the Earth without distortion, so we use many overlapping maps, each accurate locally. Smooth manifolds generalize this idea: complex shapes modeled locally by simple, flat coordinates. Concept: A smooth manifold is a space that, around every point, admits a coordinate chart mapping a neighborhood to an open set of Euclidean space R^n. An atlas is a collection of such charts whose overlaps are connected by smooth transition maps. This smooth compatibility lets us differentiate functions and define geometric objects (like tangent vectors) consistently, independent of which chart we use. Example: The 2-sphere S^2 has no single distortion-free map. Two stereographic charts (from the north and south poles) cover S^2. On their overlap, the transition map is smooth; therefore S^2 is a smooth 2-manifold. Tangent spaces at points on S^2 are 2D planes touching the sphere, and we can compute tangent vectors either from derivatives of the stereographic coordinate maps or as vectors orthogonal to the radius vector at the point.

02Intuition & Analogies

Hook: Think of hiking. Standing anywhere on a mountain, your immediate surroundings feel flat—you can walk north, south, east, or west without curving noticeably. But globally, the mountain is curved. This “locally flat, globally curved” idea captures manifolds. Concept: A chart is like a local hiking map: it gives you x–y coordinates for points near where you stand. Different charts (maps) cover different parts of the mountain, and where they overlap you can switch coordinates according to a transition rule. If those switching rules are smooth (no sharp corners or tears), then you can take derivatives, measure velocities, and do calculus consistently. Example: Picture the Earth with two big maps, one centered on the Atlantic and one on the Pacific. Where they overlap, you can translate a location’s coordinates from one map to the other smoothly. A tangent vector is the arrow representing your instantaneous walking direction and speed. On one map its components might be (3,1) (3 units east, 1 unit north). On the other map, the same physical direction becomes new numbers computed by a smooth conversion rule (the Jacobian of the transition map). On a sphere in 3D, you can also notice that “tangent arrows” are exactly the ones perpendicular to the radius from the center—arrows that don’t make you leave the surface immediately.

03Formal Definition

Hook: To make calculus work on curved spaces, we formalize “local coordinates” and “smooth change of coordinates.” Concept: An n-dimensional topological manifold M is a Hausdorff, second-countable space in which every point p ∈ M has a neighborhood U and a homeomorphism (chart)

φ :

U \to V

⊂

R^{n}

. An atlas is a collection of charts { (

U_{i}

φ_{i})

} whose domains cover M. The atlas is smooth (C^∞) if for any overlap

U_{i}

∩

U_{j}

≠ ∅, the transition maps

φ_{j}

∘ φ_

i^{- 1}

: φ_i(

U_{i}

∩

U_{j}

) → φ_j(

U_{i}

∩

U_{j}

) are C^∞ diffeomorphisms onto their images. A smooth structure is a maximal smooth atlas containing all charts compatible with the given atlas. The tangent space

T_{p}

M can be defined in multiple equivalent ways: (1) as derivations at p—linear maps D:

C^{\infty} (M) \to R

satisfying Leibniz’s rule; (2) as equivalence classes of smooth curves γ with

γ (0) = p,

where

γ_{1}

∼

γ_{2}

if their coordinate representations have the same derivative at 0; or (3) via pushforward of basis vectors through a chart inverse. If F:

M \to N

is smooth, its differential (pushforward) at p is a linear map d

F_{p}

T_{p}

M →

T_{F (p)}

N, acting in coordinates via Jacobians and satisfying the chain rule.

04When to Use

Hook: Whenever you need calculus on curved spaces—optics on lenses, trajectories on curved surfaces, shape analysis—you need manifolds and tangent spaces. Concept: Use smooth manifolds to model spaces that are locally Euclidean but globally curved or constrained. Tangent spaces capture first-order behavior: velocities of curves, gradients restricted to surfaces, and linear approximations of maps. Use cases:

Geometry and physics: Describing motion on a surface, general relativity (spacetime as a 4D manifold), and field theories.
Machine learning: Optimization on constraint sets (e.g., unit sphere, Stiefel/Grassmann manifolds), manifold learning, and normalizing flows on manifolds.
Computer graphics and robotics: Surface parametrizations, texture mapping (charts), and computing directions along surfaces (tangent frames).
Numerical simulation: Constrained differential equations (DAEs) solved by projecting velocities onto tangent spaces; coordinate changes to avoid singularities.
Software design: Implement charts and Jacobians to transform vectors and covectors between coordinate systems robustly, using transition maps for consistency.

⚠️Common Mistakes

Hook: Most confusion comes from mixing up points, coordinates, and vectors—and forgetting that vectors depend on the chart. Concept and how to avoid errors:

Confusing a point with its coordinates: A point p ∈ M is not the same as φ(p) ∈ R^n. Always label which chart you’re using and convert carefully.
Ignoring transition maps: Components of a tangent vector change with coordinates via a Jacobian. If you add vectors from different charts without converting, results are meaningless.
Assuming a single chart covers all of M: Many manifolds (like S^2) need multiple charts. Use atlases to avoid singularities (e.g., stereographic charts exclude one pole each).
Mixing ambient and intrinsic definitions: On embedded manifolds, it’s tempting to treat any R^m vector as tangent. Check tangency (e.g., v ⋅ n = 0 for spheres) or compute using chart Jacobians.
Numerical pitfalls: Finite-difference Jacobians can be unstable near singular points (e.g., at the excluded pole). Use small but not tiny step sizes, normalize where needed, and test overlap regions away from singularities.
Forgetting smoothness: Homeomorphisms are not enough; require C^∞-smooth transition maps to claim a smooth structure.

Key Formulas

Local Euclidean property

Manifold: \forall p \in M, \exists (U, φ), p \in U, φ : U \to V \subset R^{n} homeomorphism .

Explanation: Every point has a neighborhood that looks like open Euclidean space via a chart. This enables doing local calculus on the manifold.

Smooth transition maps

Compatibility: φ_{j} \circ φ_{i}^{- 1} : φ_{i} (U_{i} \cap U_{j}) \to φ_{j} (U_{i} \cap U_{j}) is C^{\infty} .

Explanation: Charts are smoothly compatible if their overlap maps are smooth. This guarantees derivatives are consistent across charts.

Derivation definition of tangent space

T_{p} M = {D : C^{\infty} (M) \to R ∣ D linear, D (f g) = f (p) D (g) + g (p) D (f)} .

Explanation: A tangent vector is a derivation at p. This abstract view does not depend on an embedding and works on any smooth manifold.

Curve definition (equivalence classes)

T_{p} M ≅ {[γ] ∣ γ : (- ϵ, ϵ) \to M, γ (0) = p}

Explanation: Two curves are equivalent if their coordinate derivatives at 0 match under any (hence every) chart. Their equivalence class gives a tangent vector.

Pushforward

d F_{p} : T_{p} M \to T_{F (p)} N, [γ]^{'} (0) \mapsto [F \circ γ]^{'} (0) .

Explanation: The differential of F sends velocities along curves through F. In coordinates it is multiplication by the Jacobian of F’s coordinate representation.

Chain rule

D (g \circ f) (x) = D g (f (x)) D f (x) .

Explanation: The Jacobian of a composition is the product of Jacobians. This underlies how components change between charts and how pushforwards compose.

Coordinate change for tangent vectors

[v]_{ψ} = D (ψ \circ φ^{- 1}) (φ (p)) [v]_{φ} .

Explanation: Given vector components in chart $φ,$ you obtain components in chart ψ by multiplying by the Jacobian of the transition map at p.

Inverse stereographic (north)

σ_{N}^{- 1} (u, v) = (\frac{2 u}{u ^{2} + v ^{2} + 1}, \frac{2 v}{u ^{2} + v ^{2} + 1}, \frac{u ^{2} + v ^{2} - 1}{u ^{2} + v ^{2} + 1}) .

Explanation: This maps plane coordinates to points on $S^{2}$ excluding the north pole. It provides a concrete chart inverse for computations.

Stereographic transition map

τ (u, v) = (\frac{u}{u ^{2} + v ^{2}}, \frac{v}{u ^{2} + v ^{2}})

Explanation: The transition from the north to south stereographic chart is inversion through the unit circle. It is smooth away from the origin (excluded pole).

Tangent plane to the sphere

T_{p} S^{2} = {v \in R^{3} ∣ v \cdot p = 0} .

Explanation: At a point p on the unit sphere, tangent vectors are exactly those orthogonal to p. This links intrinsic and embedded definitions.

Dimension of tangent space

dim T_{p} M = n

Explanation: The tangent space at any point of an n-manifold has n degrees of freedom. Locally, it behaves like $R^{n}$ .

Complexity Analysis

The computational tasks here rely primarily on evaluating Jacobians via finite differences and performing small linear algebra operations. For a function f:

R^{n}

→

R^{m}

, a central-difference Jacobian requires 2 function evaluations per input dimension, for a total of 2n evaluations of an m-vector. If the cost of a single evaluation is

C_{f}

, the total time is O(n ·

C_{f}

+ m n) to assemble the matrix (the m n term accounts for bookkeeping). Memory usage is O(m n) to store the Jacobian, plus O(n) for temporary vectors. In the stereographic examples, maps are closed-form and cheap (

C_{f}

≈ O(1)), so the Jacobian computation is dominated by O(m n) arithmetic. Transition maps between

R^{2}

charts produce 2×2 Jacobians; arithmetic is constant time, and computing determinants or inverses remains O(1). For embedded tangent computations (J:

R^{2}

→

R^{3}

), forming J is O(1) per evaluation, and projecting a 3D vector into chart coordinates via the normal equations involves forming

G = J^{T}

J (2×2), inverting G (O(1)), and multiplying matrices/vectors (O(1)). Numerical stability considerations dominate practical performance. Very small finite-difference steps h can cause catastrophic cancellation; too large h reduces accuracy through truncation error. A typical heuristic is h ≈ 1e−6 to 1e−5 for double precision, away from singularities. Around chart singularities (e.g., near excluded poles for stereographic projection), both function values and Jacobians can become ill-conditioned; complexity remains the same but accuracy degrades, so code should avoid those regions or switch charts. Overall, these routines run in constant time for fixed small dimensions and scale linearly with n and m for general

R^{n}

→

R^{m}

functions.

Code Examples

Numeric Jacobian of a Transition Map (Stereographic North → South)

1 #include <bits/stdc++.h>
2 using namespace std;
3 
4 using Vec = vector<double>;
5 using Mat = vector<vector<double>>;
6 
7 // Central-difference Jacobian for f: R^n -> R^m
8 Mat jacobian(const function<Vec(const Vec&)>& f, const Vec& x, double h = 1e-6) {
9     size_t n = x.size();
10     Vec f0 = f(x);
11     size_t m = f0.size();
12     Mat J(m, vector<double>(n, 0.0));
13     for (size_t j = 0; j < n; ++j) {
14         Vec xp = x, xm = x;
15         xp[j] += h; xm[j] -= h;
16         Vec fp = f(xp);
17         Vec fm = f(xm);
18         for (size_t i = 0; i < m; ++i) {
19             J[i][j] = (fp[i] - fm[i]) / (2.0 * h);
20         }
21     }
22     return J;
23 }
24 
25 // Transition map tau: R^2 -> R^2 between stereographic charts (north -> south)
26 // tau(u,v) = (u/(u^2+v^2), v/(u^2+v^2))
27 Vec tau(const Vec& uv) {
28     double u = uv[0], v = uv[1];
29     double r2 = u*u + v*v;
30     // Origin (0,0) is excluded (corresponds to the removed pole).
31     return { u / r2, v / r2 };
32 }
33 
34 int main() {
35     cout.setf(std::ios::fixed); cout << setprecision(6);
36 
37     // Point in the overlap (avoid the origin):
38     Vec uv = {1.2, -0.7};
39 
40     // Compute Jacobian of the transition map at uv
41     Mat J = jacobian(tau, uv);
42 
43     // Print Jacobian and its determinant to check local invertibility
44     cout << "J = [" << J[0][0] << ", " << J[0][1] << "; "
45          << J[1][0] << ", " << J[1][1] << "]\n";
46     double det = J[0][0]*J[1][1] - J[0][1]*J[1][0];
47     cout << "det(J) = " << det << "\n";
48 
49     // Transform a tangent vector's components via the Jacobian
50     Vec v_phi = {0.3, 1.1}; // components in north chart
51     Vec v_psi = { J[0][0]*v_phi[0] + J[0][1]*v_phi[1],
52                   J[1][0]*v_phi[0] + J[1][1]*v_phi[1] };
53     cout << "v in south coords = (" << v_psi[0] << ", " << v_psi[1] << ")\n";
54     return 0;
55 }
56

We implement a general central-difference Jacobian and apply it to the stereographic transition map τ between the north and south charts on S^2. We then use the Jacobian to convert tangent vector components from the north-chart coordinates to the south-chart coordinates. The determinant indicates local invertibility (away from the excluded origin).

Time: O(m n) function values, here constant O(1) for 2×2 with closed-form τ.Space: O(m n) to store the Jacobian, i.e., constant O(1) here.

Tangent Basis on S^2 from a Chart Inverse (Embedded in R^3)

1 #include <bits/stdc++.h>
2 using namespace std;
3 
4 using Vec = vector<double>;
5 using Mat = vector<vector<double>>;
6 
7 Mat jacobian(const function<Vec(const Vec&)>& f, const Vec& x, double h = 1e-6) {
8     size_t n = x.size();
9     Vec f0 = f(x);
10     size_t m = f0.size();
11     Mat J(m, vector<double>(n, 0.0));
12     for (size_t j = 0; j < n; ++j) {
13         Vec xp = x, xm = x;
14         xp[j] += h; xm[j] -= h;
15         Vec fp = f(xp);
16         Vec fm = f(xm);
17         for (size_t i = 0; i < m; ++i) {
18             J[i][j] = (fp[i] - fm[i]) / (2.0 * h);
19         }
20     }
21     return J;
22 }
23 
24 // Inverse stereographic projection from north chart: (u,v) -> (x,y,z) on S^2
25 Vec phi_inv_N(const Vec& uv) {
26     double u = uv[0], v = uv[1];
27     double r2 = u*u + v*v;
28     double denom = r2 + 1.0;
29     double x = 2.0*u/denom;
30     double y = 2.0*v/denom;
31     double z = (r2 - 1.0)/denom;
32     return {x, y, z};
33 }
34 
35 // Dot product and utilities
36 double dot3(const Vec& a, const Vec& b) { return a[0]*b[0] + a[1]*b[1] + a[2]*b[2]; }
37 
38 int main() {
39     cout.setf(std::ios::fixed); cout << setprecision(6);
40 
41     // Choose coordinates away from singularities
42     Vec uv = {0.8, -0.4};
43 
44     // Point on the sphere
45     Vec p = phi_inv_N(uv);
46 
47     // Jacobian J: R^2 -> R^3 (columns are embedded tangent basis vectors)
48     Mat J = jacobian(phi_inv_N, uv);
49 
50     // Extract tangent basis vectors from columns of J
51     Vec e1 = {J[0][0], J[1][0], J[2][0]};
52     Vec e2 = {J[0][1], J[1][1], J[2][1]};
53 
54     // Verify tangency: e_i · p ≈ 0
55     cout << "p = (" << p[0] << ", " << p[1] << ", " << p[2] << ")\n";
56     cout << "e1·p ≈ " << dot3(e1, p) << ", e2·p ≈ " << dot3(e2, p) << "\n";
57 
58     // Map a coordinate vector (a,b) to an ambient tangent vector V = J * [a;b]
59     Vec v_coords = {1.0, 0.5};
60     Vec V = { e1[0]*v_coords[0] + e2[0]*v_coords[1],
61               e1[1]*v_coords[0] + e2[1]*v_coords[1],
62               e1[2]*v_coords[0] + e2[2]*v_coords[1] };
63     cout << "Ambient tangent V = (" << V[0] << ", " << V[1] << ", " << V[2] << ")\n";
64     cout << "V·p ≈ " << dot3(V, p) << " (should be ~0)\n";
65 
66     return 0;
67 }
68

We build the embedded tangent basis on S^2 from the Jacobian of the chart inverse φ_N^{-1}: R^2 → S^2 ⊂ R^3. The Jacobian’s columns are the coordinate basis vectors in R^3. We verify they are orthogonal to the position vector p (tangent-plane condition) and use J to map coordinate components to an ambient tangent vector.

Time: O(1) for fixed 3×2 Jacobian via finite differences (constant-time evaluation).Space: O(1) to store the 3×2 Jacobian and temporary vectors.

Changing Tangent Components Between Charts and Verification via Ambient Geometry

1 #include <bits/stdc++.h>
2 using namespace std;
3 
4 using Vec = vector<double>;
5 using Mat = vector<vector<double>>;
6 
7 Mat jacobian(const function<Vec(const Vec&)>& f, const Vec& x, double h = 1e-6) {
8     size_t n = x.size();
9     Vec f0 = f(x);
10     size_t m = f0.size();
11     Mat J(m, vector<double>(n, 0.0));
12     for (size_t j = 0; j < n; ++j) {
13         Vec xp = x, xm = x;
14         xp[j] += h; xm[j] -= h;
15         Vec fp = f(xp);
16         Vec fm = f(xm);
17         for (size_t i = 0; i < m; ++i) {
18             J[i][j] = (fp[i] - fm[i]) / (2.0 * h);
19         }
20     }
21     return J;
22 }
23 
24 // S^2: inverse stereographic (north) and forward stereographic (north/south)
25 Vec phi_inv_N(const Vec& uv) {
26     double u = uv[0], v = uv[1];
27     double r2 = u*u + v*v;
28     double denom = r2 + 1.0;
29     return { 2.0*u/denom, 2.0*v/denom, (r2 - 1.0)/denom };
30 }
31 
32 Vec sigma_N(const Vec& xyz) { // (x,y,z) -> (u,v) from north chart
33     double x = xyz[0], y = xyz[1], z = xyz[2];
34     double denom = 1.0 - z; // exclude z=1
35     return { x/denom, y/denom };
36 }
37 
38 Vec sigma_S(const Vec& xyz) { // (x,y,z) -> (u,v) from south chart
39     double x = xyz[0], y = xyz[1], z = xyz[2];
40     double denom = 1.0 + z; // exclude z=-1
41     return { x/denom, y/denom };
42 }
43 
44 // Transition map tau = sigma_S o phi_inv_N
45 Vec tau(const Vec& uv) {
46     return sigma_S(phi_inv_N(uv));
47 }
48 
49 // Simple 2x2 inverse for symmetric positive definite matrix
50 bool invert2x2(const Mat& A, Mat& Ainv) {
51     double a = A[0][0], b = A[0][1], c = A[1][0], d = A[1][1];
52     double det = a*d - b*c;
53     if (fabs(det) < 1e-12) return false;
54     Ainv = {{ d/det, -b/det }, { -c/det, a/det }};
55     return true;
56 }
57 
58 int main() {
59     cout.setf(std::ios::fixed); cout << setprecision(6);
60     // Pick a point in the overlap and a tangent vector in north coordinates
61     Vec uv_N = {0.6, 0.9}; // north-chart coords (avoid origin and huge radii)
62     Vec vN = {1.2, -0.4};  // components in north chart
63 
64     // 1) Convert components using the transition map Jacobian
65     Mat Jtau = jacobian(tau, uv_N);
66     Vec vS_via_tau = { Jtau[0][0]*vN[0] + Jtau[0][1]*vN[1],
67                        Jtau[1][0]*vN[0] + Jtau[1][1]*vN[1] };
68 
69     // 2) Independently verify via ambient mapping and least-squares projection
70     //    Map point to S^2, then to south coords
71     Vec p = phi_inv_N(uv_N);
72     Vec uv_S = sigma_S(p);
73 
74     // Jacobians of chart inverses at north and of south inverse (need J_S for projection)
75     // For south, we can reuse inverse stereographic but with z sign flipped. Instead, get J of psi^{-1} via composing a small lambda.
76     auto phi_inv_S = [](const Vec& uv)->Vec {
77         double u = uv[0], v = uv[1];
78         double r2 = u*u + v*v; double denom = r2 + 1.0;
79         return { 2.0*u/denom, 2.0*v/denom, -(r2 - 1.0)/denom };
80     };
81 
82     Mat JN = jacobian(phi_inv_N, uv_N); // 3x2
83     Mat JS = jacobian(function<Vec(const Vec&)>(phi_inv_S), uv_S); // 3x2
84 
85     // Ambient tangent vector: V = JN * vN
86     Vec V = { JN[0][0]*vN[0] + JN[0][1]*vN[1],
87               JN[1][0]*vN[0] + JN[1][1]*vN[1],
88               JN[2][0]*vN[0] + JN[2][1]*vN[1] };
89 
90     // Recover south components by solving min_w ||JS*w - V||^2 via normal equations
91     Mat G = { { JS[0][0]*JS[0][0] + JS[1][0]*JS[1][0] + JS[2][0]*JS[2][0],
92                 JS[0][0]*JS[0][1] + JS[1][0]*JS[1][1] + JS[2][0]*JS[2][1] },
93               { 0.0, 0.0 } };
94     G[1][0] = G[0][1];
95     G[1][1] = JS[0][1]*JS[0][1] + JS[1][1]*JS[1][1] + JS[2][1]*JS[2][1];
96 
97     Vec b = { JS[0][0]*V[0] + JS[1][0]*V[1] + JS[2][0]*V[2],
98               JS[0][1]*V[0] + JS[1][1]*V[1] + JS[2][1]*V[2] };
99 
100     Mat Ginv;
101     bool ok = invert2x2(G, Ginv);
102     if (!ok) { cerr << "Ill-conditioned system (near singular chart).\n"; return 0; }
103     Vec vS_via_proj = { Ginv[0][0]*b[0] + Ginv[0][1]*b[1],
104                         Ginv[1][0]*b[0] + Ginv[1][1]*b[1] };
105 
106     cout << "vS via tau    = (" << vS_via_tau[0] << ", " << vS_via_tau[1] << ")\n";
107     cout << "vS via proj   = (" << vS_via_proj[0] << ", " << vS_via_proj[1] << ")\n";
108     cout << "difference    = (" << (vS_via_tau[0]-vS_via_proj[0]) << ", "
109          << (vS_via_tau[1]-vS_via_proj[1]) << ")\n";
110 
111     return 0;
112 }
113

We convert tangent vector components from the north chart to the south chart in two ways: (1) by multiplying with the Jacobian of the transition map τ, and (2) by mapping the vector into the ambient space using J_N, then solving a 2×2 least-squares system to recover components in the south chart using J_S. The two answers agree (up to numerical error), verifying the coordinate-change rule and the embedded interpretation of tangent vectors.

Time: O(1) for fixed dimensions: a few Jacobian evaluations (3×2 and 3×2), small matrix multiplications, and a 2×2 inverse.Space: O(1) for small matrices and vectors.

1	#include <bits/stdc++.h>
2	using namespace std;
3
4	using Vec = vector<double>;
5	using Mat = vector<vector<double>>;
6
7	// Central-difference Jacobian for f: R^n -> R^m
8	Mat jacobian(const function<Vec(const Vec&)>& f, const Vec& x, double h = 1e-6) {
9	size_t n = x.size();
10	Vec f0 = f(x);
11	size_t m = f0.size();
12	Mat J(m, vector<double>(n, 0.0));
13	for (size_t j = 0; j < n; ++j) {
14	Vec xp = x, xm = x;
15	xp[j] += h; xm[j] -= h;
16	Vec fp = f(xp);
17	Vec fm = f(xm);
18	for (size_t i = 0; i < m; ++i) {
19	J[i][j] = (fp[i] - fm[i]) / (2.0 * h);
20	}
21	}
22	return J;
23	}
24
25	// Transition map tau: R^2 -> R^2 between stereographic charts (north -> south)
26	// tau(u,v) = (u/(u^2+v^2), v/(u^2+v^2))
27	Vec tau(const Vec& uv) {
28	double u = uv[0], v = uv[1];
29	double r2 = uu + vv;
30	// Origin (0,0) is excluded (corresponds to the removed pole).
31	return { u / r2, v / r2 };
32	}
33
34	int main() {
35	cout.setf(std::ios::fixed); cout << setprecision(6);
36
37	// Point in the overlap (avoid the origin):
38	Vec uv = {1.2, -0.7};
39
40	// Compute Jacobian of the transition map at uv
41	Mat J = jacobian(tau, uv);
42
43	// Print Jacobian and its determinant to check local invertibility
44	cout << "J = [" << J[0][0] << ", " << J[0][1] << "; "
45	<< J[1][0] << ", " << J[1][1] << "]\n";
46	double det = J[0][0]J[1][1] - J[0][1]J[1][0];
47	cout << "det(J) = " << det << "\n";
48
49	// Transform a tangent vector's components via the Jacobian
50	Vec v_phi = {0.3, 1.1}; // components in north chart
51	Vec v_psi = { J[0][0]v_phi[0] + J[0][1]v_phi[1],
52	J[1][0]v_phi[0] + J[1][1]v_phi[1] };
53	cout << "v in south coords = (" << v_psi[0] << ", " << v_psi[1] << ")\n";
54	return 0;
55	}
56