Stereographic projection

Joshua R. Davis, jdavis@carleton.edu, revised 2020 October 2

Two-dimensional plane and sphere

First we lay out some infrastructure. Don’t worry about the details of how it works. The cart function converts polar coordinates to Cartesian coordinates. The param function is a parametrization of the punctured sphere; its inverse is the stereographic projection.

Clear[cart, param];
cart[{r_, theta_}] := 
{r Cos[theta], r Sin[theta]};
param[{u_,v_}]:=
{(2 u)/(u^2+v^2+1),(2 v)/(u^2+v^2+1),(u^2+v^2-1)/(u^2+v^2+1)};

Here’s an illustration of what the parametrization does. It wraps the entire plane around the sphere, leaving only the north pole uncovered. On a computer, we can’t show the entire plane, so we restrict ourselves to a disk.

a=5;
polarPlanePlot = ParametricPlot3D[Append[cart[{r,theta}], 0],{r,0,5},{theta,0,2 Pi}, PlotRange->{{-a, a}, {-a, a}, {-1, 1}}];
polarSpherePlot = ParametricPlot3D[param[cart[{r, theta}]],{r,0,5},{theta,0,2 Pi}, PlotRange->{{-1, 1}, {-1, 1}, {-1, 1}}];
polarPlanePlot -> polarSpherePlot

Here’s a manipulable version, which intuitively animates how the disk gets wrapped around the sphere. (Technically, the animation is a homotopy between (when u = 0) the identity function from the disk to (R^3) and (when u = 1) the param map from the disk to (R^3).)

a=5;
Manipulate[
ParametricPlot3D[
(1 - u) Append[cart[{r,theta}], 0] + u param[cart[{r, theta}]],
{r,0,5},{theta,0,2 Pi}, PlotRange->{{-a, a}, {-a, a}, {-1, 1}}], 
{u, 0, 1}]

Here’s a different manipulable, where you can vary the radius of the disk and see how the coverage of the sphere changes. Notice that a disk of radius 1 covers exactly a hemisphere.

a = 5;
Manipulate[
Module[{polarPlanePlot, polarSpherePlot}, polarPlanePlot = ParametricPlot3D[Append[cart[{r,theta}], 0],{r,0,r0},{theta,0,2 Pi}, PlotRange->{{-a, a}, {-a, a}, {-1, 1}}];
polarSpherePlot = ParametricPlot3D[param[cart[{r, theta}]],{r,0,r0},{theta,0,2 Pi}, PlotRange->{{-1, 1}, {-1, 1}, {-1, 1}}];
polarPlanePlot -> polarSpherePlot], 
{r0, 1, 5}]

Here’s the same idea, but using a rectilinear grid instead of a polar grid. The image of the grid on the sphere is pretty, I think.

a = 5;
cartPlanePlot = ParametricPlot3D[{u,v, 0},{u,-a,a},{v,-a,a}, PlotRange->{{-a, a}, {-a, a}, {-1, 1}}];
cartSpherePlot = ParametricPlot3D[param[{u,v}],{u,-a,a},{v,-a,a}, PlotRange->{{-1, 1}, {-1, 1}, {-1, 1}}];
cartPlanePlot -> cartSpherePlot

Notice that the grid lines are perpendicular on the sphere, just as they are perpendicular on the plane. Also notice that the grid squares all have the same area on the plane, but they have wildly varying areas on the sphere.

Checking inverses in arbitrary dimension

Here’s some infrastructure for making stereographic projections in arbitrary dimension. Don’t worry about the details of how it works.

Clear[fExpr, fInvExpr, x, y];
fExpr[n_] := 
Module[{xx, most}, 
xx = Array[x, n];
most = Table[2 x[k] / (xx.xx + 1), {k, 1, n}];
Append[most, (xx.xx - 1) / (xx.xx + 1)]];
fInvExpr[n_] := 
Table[y[k] / (1 - y[n + 1]), {k, 1, n}];

These are the expressions for f (the parametrization) and (f^(-1)) (the projection) used above.

fExpr[2]
fInvExpr[2]

Here’s the n = 3 case.

n = 3;
f = fExpr[n]
fInv = fInvExpr[n]

Intuitively, f outputs y-values and (f^(-1)) inputs y-values. So we should be able to plug the outputs of f into the inputs of (f^(-1)) and end up with the identity function id(x) = x. Yep.

fInvF = fInv /. Table[y[k] -> f[[k]], {k, 1, n + 1}]
Simplify[fInvF]

Similarly, (f^(-1)) outputs x-values and f inputs x-values, so we should get the identity function id(y) = y when we plug the outputs of (f^(-1)) into the inputs of f. Do we?

fFInv = f /. Table[x[k] -> fInv[[k]], {k, 1, n}]
Simplify[fFInv]

The problem in that last computation is that Mathematica doesn’t know a crucial piece of information: The only y that we’re considering are y of norm 1. We can either simplify manually ourselves, or we can give that extra information to Mathematica’s Simplify function. Either way, we get id(y) = y.

normSqY = Table[y[k], {k, 1, n + 1}].Table[y[k], {k, 1, n + 1}]
Simplify[fFInv, {normSqY == 1}]

Try re-running all this code starting with a different choice of n. In all cases, you should be reassured that f and (f^(-1)) really are inverses.

Jacobian

Here is some Mathematica code for computing the Jacobian matrix of a function f : (R^n)-> (R^m). It has nothing to do with stereographic projection.

Clear[jac];
jac[f_, vars_] := 
Transpose[Table[D[f, vars[[k]]], {k, 1, Length[vars]}]];

Here’s a demonstration.

f = {Cos[x y], y, x^3 Sin[y]}
jacF = jac[f, {x, y}]
MatrixForm[jacF]

Here’s the Jacobian when f is the parametrization of the 2-sphere by the inverse of stereographic projection.

f = fExpr[2]
jacF = Simplify[jac[f, {x[1], x[2]}]]
MatrixForm[jacF]

Here’s the transpose Jacobian times the Jacobian. Have you ever seen this matrix in this course?

jacFTJacF = Simplify[Transpose[jacF].jacF]
MatrixForm[jacFTJacF]