2336-Finding area preserving transformations.sagews

Suppose we'd like to randomly generate some points that are uniformly distributed throughout the unit disk.  Since the disk is known to be parametrized by

$$

\begin{array}{lll}

x(r,\theta) & = & r \cos(\theta) \\

y(r,\theta) & = & r \sin(\theta),

\end{array}

$$

one natural guess might be to pick $r$ uniformly from the unit interval and $\theta$ uniformly from $[0,2\pi)$.  We could then determine $x$ and $y$ from the above equations.  We can code this as follows.

def pt_in_disk(): 
    r = random()
    t = 2*pi*random()
    return [n(r*cos(t)), n(r*sin(t))]

Let's try it.

pts = [pt_in_disk() for i in xrange(10000)]
show(point(pts) + circle((0,0),1, thickness=4),
    aspect_ratio=1)

Uh-oh - the points are clearly not uniformly distributed; they appear to cluster too tightly near the center.  In retrospect, the problem is obvious: The circle with radius 1/2 contains just about half the points but only a quarter of the area.

len(filter(lambda v: norm(v)<1/2, [vector(pt) for pt in pts]))

5018

At MathWorld (http://mathworld.wolfram.com/DiskPointPicking.html) we see that we should pick $r$ uniformly from the unit interval and $\theta$ uniformly from $[0,2\pi)$ and then determine $x$ and $y$ according the equations

$$

\begin{array}{lll}

x(r,\theta) & = & \sqrt{r} \cos(\theta) \\

y(r,\theta) & = & \sqrt{r} \sin(\theta).

\end{array}

$$

Let's try that.

def pt_in_disk(): 
    r = random()
    t = 2*pi*random()
    return [n(sqrt(r)*cos(t)), n(sqrt(r)*sin(t))]
pts = [pt_in_disk() for i in xrange(10000)]
show(point(pts) + circle((0,0),1, thickness=4),
    aspect_ratio=1)

Much better. Ultimately, this works because the determinant of the Jacobian of the transformation $T(r,\theta )=\sqrt{r}(\cos (\theta ),\sin (\theta ))$ is constant.

r,t = var('r,t')
T = [sqrt(r)*cos(t), sqrt(r)*sin(t)]
jac = det(jacobian(T,[r,t]))
jac.trig_simplify()

1/2

We can use this criterion on other transformations.  For example, this parametrization of an ellipse has a Jacobian with constant determinant.

a,b = var('a,b')
T = [a*sqrt(r)*cos(t), b*sqrt(r)*sin(t)]
jac = det(jacobian(T,[r,t]))
jac.trig_simplify()

1/2*a*b

Thus, we can use that transformation to uniformly generate points in an ellipse.

def pt_in_disk(): 
    r = random()
    t = 2*pi*random()
    return [n(8*sqrt(r)*cos(t)), n(5*sqrt(r)*sin(t))]
pts = [pt_in_disk() for i in xrange(10000)]
show(point(pts) + parametric_plot([8*cos(t),5*sin(t)],[t,0,2*pi], color='black', thickness=4),
    aspect_ratio=1)

Looks good!

Finding area preserving transformations

Of course, the big question is, how might we find such transformations?  One technique is to start with a known parametrization and scale the variables by some non-linear function.  Typically, that function will need to solve some differential equation to work.  Suppose, for example, we start with the well-known parametrization of a disk.

$$

\begin{array}{lll}

x(r,\theta) & = & r \cos(\theta) \\

y(r,\theta) & = & r \sin(\theta),

\end{array}

$$

We suspect that some radially re-scaled transformation might have constant Jacobian determinant.  Let's write that down in terms of some, as yet, unknown function.

$$

\begin{array}{lll}

x(r,\theta) & = & f(r) \cos(\theta) \\

y(r,\theta) & = & f(r) \sin(\theta),

\end{array}

$$

We now compute the determinant of the Jacobian of this transformation.

r,t = var('r,t')
f = function('f',r)
T = [f*cos(t), f*sin(t)]
jac = det(jacobian(T,[r,t]))
jac.trig_simplify()

f(r)*D[0](f)(r)

If we want this to be constant, then we must solve the differential equation $f(r)f'(r)=c$.  Any constant will do, as long as the differential equation can be solved, but note that $f(r)=\sqrt{r}$ solves the equation with $c=1/2$.

Uniform points in the sphere

Let's now consider the uniform distribution of points throughout a sphere, which is quite a bit harder.  The basic transformation is given by the spherical coordinate formulae.  We can generalize to the case of an ellipsoid by simply multiplying the $x$, $y$, and $z$ equations by $a$, $b$, and $c$.

a,b,c, rho,phi,theta = var('a,b,c, rho,phi,theta')
T = [a*rho*sin(phi)*cos(theta), b*rho*sin(phi)*sin(theta), c*rho*cos(phi)]
jac = det(jacobian(T,[rho,phi,theta]))
jac.trig_simplify()

a*b*c*rho^2*sin(phi)

The $\rho^2 \sin(\phi)$ should be familiar to anyone who's worked with cylindrical coordinates.  Let's now try scaling both the $\rho$ and $\phi$ variables with unknown functions $f$ and $g$.

a,b,c, rho,phi,theta = var('a,b,c, rho,phi,theta')
f = function('f', rho)
g = function('g', phi)
T = [a*rho*sin(phi)*cos(theta), b*rho*sin(phi)*sin(theta), c*rho*cos(phi)]
T = [expr(rho=f, phi=g) for expr in T]
jac = det(jacobian(T, [rho,phi,theta]))
jac.trig_simplify()

a*b*c*sin(g(phi))*f(rho)^2*D[0](f)(rho)*D[0](g)(phi)

Looking at this, I would say that we need to solve the differential equations $f(\rho)^2 f'(\rho) = c_1$ and $\sin(g(\phi)) g'(\phi) = c_2$.  Furthermore, we'll need $f(0)=g(0)=0$.  The first equation is easy.

desolve(f^2 * diff(f,rho)==1,f, ics=[0,0])

1/3*f(rho)^3 == rho

The second equation has a couple of complications.  First, experimentation led me to realize that the choice of the constant $c_2$ is a bit important down the line.  I ended up choosing $c_2 = 1/\pi$.

desolve(sin(g)*diff(g,phi) == 1/pi, g, ics = [0,0])

-pi*cos(g(phi)) == -pi + phi

So I guess we'll take $g(\phi) = \arccos((\phi-\pi)/\pi) = \sqrt{1-{(\phi -\pi )^2}/{\pi ^2}}$.  Note that this is never negative; presumably, we lost some information applying the arccosine function.  We can still try it and run a numerical check at the end.  First, let's check the determinant of the Jacobian to see if it is, indeed, constant.

T = [a*rho*sin(phi)*cos(theta), b*rho*sin(phi)*sin(theta), c*rho*cos(phi)]
T = [expr(phi=arccos((pi-phi)/pi), rho=rho^(1/3)) for expr in T]
jac = det(jacobian(T, [rho,phi,theta]))
jac.trig_simplify()

1/3*a*b*c/pi

OK, let's try it on the sphere.  Note the (2*randint(0,1)-1) construct,which slaps a $\pm 1$ in front of the $z$-coordinate.

a=1; b=1; c=1;
def pt_in_sphere(): 
    rho = random()
    phi = pi*random()
    theta = 2*pi*random()
    return [sqrt(-(pi - phi)^2/pi^2 + 1)*a*rho^(1/3)*cos(theta), \
    sqrt(-(pi - phi)^2/pi^2 + 1)*b*rho^(1/3)*sin(theta), (2*randint(0,1)-1)*(pi - phi)*c*rho^(1/3)/pi]
pts = [pt_in_sphere() for i in xrange(10000)]

The following computations should yield something close to $1250$, since approximately $1/8$ of the points should be inside a sub-sphere of radius $1/2$.

len(filter(lambda v: norm(v)<1/2, [vector(pt)-vector([0.1,0.2,0.1]) for pt in pts]))

1251

You can try to visualize, too, but it might look silly.

point(pts)