Riemann Mapping
AUTHORS:
Development supported by NSF award No. 0702939.
Bases: object
The Riemann_Map class computes an interior or exterior Riemann map, or an Ahlfors map of a region given by the supplied boundary curve(s) and center point. The class also provides various methods to evaluate, visualize, or extract data from the map.
A Riemann map conformally maps a simply connected region in the complex plane to the unit disc. The Ahlfors map does the same thing for multiply connected regions.
Note that all the methods are numerical. As a result all answers have some imprecision. Moreover, maps computed with small number of collocation points, or for unusually shaped regions, may be very inaccurate. Error computations for the ellipse can be found in the documentation for analytic_boundary() and analytic_interior().
[BSV] provides an overview of the Riemann map and discusses the research that lead to the creation of this module.
INPUT:
The following inputs may be passed in as named parameters:
The following inputs may be passed as named parameters in unusual circumstances:
EXAMPLES:
The unit circle identity map:
sage: f(t) = e^(I*t)
sage: fprime(t) = I*e^(I*t)
sage: m = Riemann_Map([f], [fprime], 0) # long time (4 sec)
sage: m.plot_colored() + m.plot_spiderweb() # long time
Graphics object consisting of 22 graphics primitives
The exterior map for the unit circle:
sage: m = Riemann_Map([f], [fprime], 0, exterior=True) # long time (4 sec)
sage: #spiderwebs are not supported for exterior maps
sage: m.plot_colored() # long time
Graphics object consisting of 1 graphics primitive
The unit circle with a small hole:
sage: f(t) = e^(I*t)
sage: fprime(t) = I*e^(I*t)
sage: hf(t) = 0.5*e^(-I*t)
sage: hfprime(t) = 0.5*-I*e^(-I*t)
sage: m = Riemann_Map([f, hf], [fprime, hfprime], 0.5 + 0.5*I)
sage: #spiderweb and color plots cannot be added for multiply
sage: #connected regions. Instead we do this.
sage: m.plot_spiderweb(withcolor = True) # long time
Graphics object consisting of 3 graphics primitives
A square:
sage: ps = polygon_spline([(-1, -1), (1, -1), (1, 1), (-1, 1)])
sage: f = lambda t: ps.value(real(t))
sage: fprime = lambda t: ps.derivative(real(t))
sage: m = Riemann_Map([f], [fprime], 0.25, ncorners=4)
sage: m.plot_colored() + m.plot_spiderweb() # long time
Graphics object consisting of 22 graphics primitives
Compute rough error for this map:
sage: x = 0.75 # long time
sage: print "error =", m.inverse_riemann_map(m.riemann_map(x)) - x # long time
error = (-0.000...+0.0016...j)
A fun, complex region for demonstration purposes:
sage: f(t) = e^(I*t)
sage: fp(t) = I*e^(I*t)
sage: ef1(t) = .2*e^(-I*t) +.4+.4*I
sage: ef1p(t) = -I*.2*e^(-I*t)
sage: ef2(t) = .2*e^(-I*t) -.4+.4*I
sage: ef2p(t) = -I*.2*e^(-I*t)
sage: pts = [(-.5,-.15-20/1000),(-.6,-.27-10/1000),(-.45,-.45),(0,-.65+10/1000),(.45,-.45),(.6,-.27-10/1000),(.5,-.15-10/1000),(0,-.43+10/1000)]
sage: pts.reverse()
sage: cs = complex_cubic_spline(pts)
sage: mf = lambda x:cs.value(x)
sage: mfprime = lambda x: cs.derivative(x)
sage: m = Riemann_Map([f,ef1,ef2,mf],[fp,ef1p,ef2p,mfprime],0,N = 400) # long time
sage: p = m.plot_colored(plot_points = 400) # long time
ALGORITHM:
This class computes the Riemann Map via the Szego kernel using an adaptation of the method described by [KT].
REFERENCES:
[KT] | N. Kerzman and M. R. Trummer. “Numerical Conformal Mapping via the Szego kernel”. Journal of Computational and Applied Mathematics, 14(1-2): 111–123, 1986. |
[BSV] | M. Bolt, S. Snoeyink, E. Van Andel. “Visual representation of the Riemann map and Ahlfors map via the Kerzman-Stein equation”. Involve 3-4 (2010), 405-420. |
Computes the Riemann map on a grid of points. Note that these points are complex of the form z = x + y*i.
INPUT:
OUTPUT:
EXAMPLES:
General usage:
sage: f(t) = e^(I*t) - 0.5*e^(-I*t)
sage: fprime(t) = I*e^(I*t) + 0.5*I*e^(-I*t)
sage: m = Riemann_Map([f], [fprime], 0)
sage: data = m.compute_on_grid([],5)
sage: print data[0][8,1]
(-0.0879...+0.9709...j)
Returns a discretized version of the Szego kernel for each boundary function.
INPUT:
The following inputs may be passed in as named parameters:
OUTPUT:
A list of points of the form [t value, value of the Szego kernel at that t].
EXAMPLES:
Generic use:
sage: f(t) = e^(I*t) - 0.5*e^(-I*t)
sage: fprime(t) = I*e^(I*t) + 0.5*I*e^(-I*t)
sage: m = Riemann_Map([f], [fprime], 0)
sage: sz = m.get_szego(boundary=0)
sage: points = m.get_szego(absolute_value=True)
sage: list_plot(points)
Graphics object consisting of 1 graphics primitive
Extending the points by a spline:
sage: s = spline(points)
sage: s(3*pi / 4)
0.0012158...
sage: plot(s,0,2*pi) # plot the kernel
Graphics object consisting of 1 graphics primitive
The unit circle with a small hole:
sage: f(t) = e^(I*t)
sage: fprime(t) = I*e^(I*t)
sage: hf(t) = 0.5*e^(-I*t)
sage: hfprime(t) = 0.5*-I*e^(-I*t)
sage: m = Riemann_Map([f, hf], [fprime, hfprime], 0.5 + 0.5*I)
Getting the szego for a specifc boundary:
sage: sz0 = m.get_szego(boundary=0)
sage: sz1 = m.get_szego(boundary=1)
Returns an array of points of the form [t value, theta in e^(I*theta)], that is, a discretized version of the theta/boundary correspondence function. In other words, a point in this array [t1, t2] represents that the boundary point given by f(t1) is mapped to a point on the boundary of the unit circle given by e^(I*t2).
For multiply connected domains, get_theta_points will list the points for each boundary in the order that they were supplied.
INPUT:
The following input must all be passed in as named parameters:
OUTPUT:
A list of points of the form [t value, theta in e^(I*theta)].
EXAMPLES:
Getting the list of points, extending it via a spline, getting the points for only the outside of a multiply connected domain:
sage: f(t) = e^(I*t) - 0.5*e^(-I*t)
sage: fprime(t) = I*e^(I*t) + 0.5*I*e^(-I*t)
sage: m = Riemann_Map([f], [fprime], 0)
sage: points = m.get_theta_points()
sage: list_plot(points)
Graphics object consisting of 1 graphics primitive
Extending the points by a spline:
sage: s = spline(points)
sage: s(3*pi / 4)
1.627660...
The unit circle with a small hole:
sage: f(t) = e^(I*t)
sage: fprime(t) = I*e^(I*t)
sage: hf(t) = 0.5*e^(-I*t)
sage: hfprime(t) = 0.5*-I*e^(-I*t)
sage: m = Riemann_Map([f, hf], [hf, hfprime], 0.5 + 0.5*I)
Getting the boundary correspondence for a specifc boundary:
sage: tp0 = m.get_theta_points(boundary=0)
sage: tp1 = m.get_theta_points(boundary=1)
Returns the inverse Riemann mapping of a point. That is, given pt on the interior of the unit disc, inverse_riemann_map() will return the point on the original region that would be Riemann mapped to pt. Note that this method does not work for multiply connected domains.
INPUT:
OUTPUT:
The point on the region that Riemann maps to the input point.
EXAMPLES:
Can work for different types of complex numbers:
sage: f(t) = e^(I*t) - 0.5*e^(-I*t)
sage: fprime(t) = I*e^(I*t) + 0.5*I*e^(-I*t)
sage: m = Riemann_Map([f], [fprime], 0)
sage: m.inverse_riemann_map(0.5 + sqrt(-0.5))
(0.406880...+0.3614702...j)
sage: m.inverse_riemann_map(0.95)
(0.486319...-4.90019052...j)
sage: m.inverse_riemann_map(0.25 - 0.3*I)
(0.1653244...-0.180936...j)
sage: import numpy as np
sage: m.inverse_riemann_map(np.complex(-0.2, 0.5))
(-0.156280...+0.321819...j)
Plots the boundaries of the region for the Riemann map. Note that this method DOES work for multiply connected domains.
INPUT:
The following inputs may be passed in as named parameters:
EXAMPLES:
General usage:
sage: f(t) = e^(I*t) - 0.5*e^(-I*t)
sage: fprime(t) = I*e^(I*t) + 0.5*I*e^(-I*t)
sage: m = Riemann_Map([f], [fprime], 0)
Default plot:
sage: m.plot_boundaries()
Graphics object consisting of 1 graphics primitive
Big blue collocation points:
sage: m.plot_boundaries(plotjoined=False, rgbcolor=[0,0,1], thickness=6)
Graphics object consisting of 1 graphics primitive
Generates a colored plot of the Riemann map. A red point on the colored plot corresponds to a red point on the unit disc.
INPUT:
The following inputs may be passed in as named parameters:
EXAMPLES:
Given a Riemann map m, general usage:
sage: f(t) = e^(I*t) - 0.5*e^(-I*t)
sage: fprime(t) = I*e^(I*t) + 0.5*I*e^(-I*t)
sage: m = Riemann_Map([f], [fprime], 0)
sage: m.plot_colored()
Graphics object consisting of 1 graphics primitive
Plot zoomed in on a specific spot:
sage: m.plot_colored(plot_range=[0,1,.25,.75])
Graphics object consisting of 1 graphics primitive
High resolution plot:
sage: m.plot_colored(plot_points=1000) # long time (29s on sage.math, 2012)
Graphics object consisting of 1 graphics primitive
To generate the unit circle map, it’s helpful to see what the colors correspond to:
sage: f(t) = e^(I*t)
sage: fprime(t) = I*e^(I*t)
sage: m = Riemann_Map([f], [fprime], 0, 1000)
sage: m.plot_colored()
Graphics object consisting of 1 graphics primitive
Generates a traditional “spiderweb plot” of the Riemann map. Shows what concentric circles and radial lines map to. The radial lines may exhibit erratic behavior near the boundary; if this occurs, decreasing linescale may mitigate the problem.
For multiply connected domains the spiderweb is by necessity generated using the forward mapping. This method is more computationally intensive. In addition, these spiderwebs cannot be added to color plots. Instead the withcolor option must be used.
In addition, spiderweb plots are not currently supported for exterior maps.
INPUT:
The following inputs may be passed in as named parameters:
EXAMPLES:
General usage:
sage: f(t) = e^(I*t) - 0.5*e^(-I*t)
sage: fprime(t) = I*e^(I*t) + 0.5*I*e^(-I*t)
sage: m = Riemann_Map([f], [fprime], 0)
Default plot:
sage: m.plot_spiderweb()
Graphics object consisting of 21 graphics primitives
Simplified plot with many discrete points:
sage: m.plot_spiderweb(spokes=4, circles=1, pts=400, linescale=0.95, plotjoined=False)
Graphics object consisting of 6 graphics primitives
Plot with thick, red lines:
sage: m.plot_spiderweb(rgbcolor=[1,0,0], thickness=3)
Graphics object consisting of 21 graphics primitives
To generate the unit circle map, it’s helpful to see what the original spiderweb looks like:
sage: f(t) = e^(I*t)
sage: fprime(t) = I*e^(I*t)
sage: m = Riemann_Map([f], [fprime], 0, 1000)
sage: m.plot_spiderweb()
Graphics object consisting of 21 graphics primitives
A multiply connected region with corners. We set min_mag higher to remove “fuzz” outside the domain:
sage: ps = polygon_spline([(-4,-2),(4,-2),(4,2),(-4,2)])
sage: z1 = lambda t: ps.value(t); z1p = lambda t: ps.derivative(t)
sage: z2(t) = -2+exp(-I*t); z2p(t) = -I*exp(-I*t)
sage: z3(t) = 2+exp(-I*t); z3p(t) = -I*exp(-I*t)
sage: m = Riemann_Map([z1,z2,z3],[z1p,z2p,z3p],0,ncorners=4) # long time
sage: p = m.plot_spiderweb(withcolor=True,plot_points=500, thickness = 2.0, min_mag=0.1) # long time
Returns the Riemann mapping of a point. That is, given pt on the interior of the mapped region, riemann_map will return the point on the unit disk that pt maps to. Note that this method only works for interior points; accuracy breaks down very close to the boundary. To get boundary corrospondance, use get_theta_points().
INPUT:
OUTPUT:
A complex number representing the point on the unit circle that the input point maps to.
EXAMPLES:
Can work for different types of complex numbers:
sage: f(t) = e^(I*t) - 0.5*e^(-I*t)
sage: fprime(t) = I*e^(I*t) + 0.5*I*e^(-I*t)
sage: m = Riemann_Map([f], [fprime], 0)
sage: m.riemann_map(0.25 + sqrt(-0.5))
(0.137514...+0.876696...j)
sage: I = CDF.gen()
sage: m.riemann_map(1.3*I)
(-1.56...e-05+0.989694...j)
sage: m.riemann_map(0.4)
(0.73324...+3.2...e-06j)
sage: import numpy as np
sage: m.riemann_map(np.complex(-3, 0.0001))
(1.405757...e-05+8.06...e-10j)
Provides an exact (for n = infinity) Riemann boundary correspondence for the ellipse with axes 1 + epsilon and 1 - epsilon. The boundary is therefore given by e^(I*t)+epsilon*e^(-I*t). It is primarily useful for testing the accuracy of the numerical Riemann_Map.
INPUT:
OUTPUT:
A theta value from 0 to 2*pi, corresponding to the point on the circle e^(I*theta)
TESTS:
Checking the accuracy of this function for different n values:
sage: from sage.calculus.riemann import analytic_boundary
sage: t100 = analytic_boundary(pi/2, 100, .3)
sage: abs(analytic_boundary(pi/2, 10, .3) - t100) < 10^-8
True
sage: abs(analytic_boundary(pi/2, 20, .3) - t100) < 10^-15
True
Using this to check the accuracy of the Riemann_Map boundary:
sage: f(t) = e^(I*t)+.3*e^(-I*t)
sage: fp(t) = I*e^(I*t)-I*.3*e^(-I*t)
sage: m = Riemann_Map([f], [fp],0,200)
sage: s = spline(m.get_theta_points())
sage: test_pt = uniform(0,2*pi)
sage: s(test_pt) - analytic_boundary(test_pt,20, .3) < 10^-4
True
Provides a nearly exact compuation of the Riemann Map of an interior point of the ellipse with axes 1 + epsilon and 1 - epsilon. It is primarily useful for testing the accuracy of the numerical Riemann Map.
INPUT:
TESTS:
Testing the accuracy of Riemann_Map:
sage: from sage.calculus.riemann import analytic_interior
sage: f(t) = e^(I*t)+.3*e^(-I*t)
sage: fp(t) = I*e^(I*t)-I*.3*e^(-I*t)
sage: m = Riemann_Map([f],[fp],0,200)
sage: abs(m.riemann_map(.5)-analytic_interior(.5, 20, .3)) < 10^-4
True
sage: m = Riemann_Map([f],[fp],0,2000)
sage: abs(m.riemann_map(.5)-analytic_interior(.5, 20, .3)) < 10^-6
True
Intermediate function for the integration in analytic_interior().
INPUT:
TESTS:
This is primarily tested implicitly by analytic_interior(). Here is a simple test:
sage: from sage.calculus.riemann import cauchy_kernel
sage: cauchy_kernel(.5,(.3, .1+.2*I, 10,'c'))
(-0.584136405997...+0.5948650858950...j)
Convert from a (Numpy) array of complex numbers to its corresponding matrix of RGB values. For internal use of plot_colored() only.
INPUT:
OUTPUT:
An \(N \times M \times 3\) floating point Numpy array X, where X[i,j] is an (r,g,b) tuple.
EXAMPLES:
sage: from sage.calculus.riemann import complex_to_rgb
sage: import numpy
sage: complex_to_rgb(numpy.array([[0, 1, 1000]], dtype = numpy.complex128))
array([[[ 1. , 1. , 1. ],
[ 1. , 0.05558355, 0.05558355],
[ 0.17301243, 0. , 0. ]]])
sage: complex_to_rgb(numpy.array([[0, 1j, 1000j]], dtype = numpy.complex128))
array([[[ 1. , 1. , 1. ],
[ 0.52779177, 1. , 0.05558355],
[ 0.08650622, 0.17301243, 0. ]]])
TESTS:
sage: complex_to_rgb([[0, 1, 10]])
Traceback (most recent call last):
...
TypeError: Argument 'z_values' has incorrect type (expected numpy.ndarray, got list)
Converts a grid of complex numbers into a matrix containing rgb data for the Riemann spiderweb plot.
INPUT:
OUTPUT:
An \(N x M x 3\) floating point Numpy array X, where X[i,j] is an (r,g,b) tuple.
EXAMPLES:
sage: from sage.calculus.riemann import complex_to_spiderweb
sage: import numpy
sage: zval = numpy.array([[0, 1, 1000],[.2+.3j,1,-.3j],[0,0,0]],dtype = numpy.complex128)
sage: deriv = numpy.array([[.1]],dtype = numpy.float64)
sage: complex_to_spiderweb(zval, deriv,deriv, 4,4,[0,0,0],1,False,0.001)
array([[[ 1., 1., 1.],
[ 1., 1., 1.],
[ 1., 1., 1.]],
[[ 1., 1., 1.],
[ 0., 0., 0.],
[ 1., 1., 1.]],
[[ 1., 1., 1.],
[ 1., 1., 1.],
[ 1., 1., 1.]]])
sage: complex_to_spiderweb(zval, deriv,deriv, 4,4,[0,0,0],1,True,0.001)
array([[[ 1. , 1. , 1. ],
[ 1. , 0.05558355, 0.05558355],
[ 0.17301243, 0. , 0. ]],
[[ 1. , 0.96804683, 0.48044583],
[ 0. , 0. , 0. ],
[ 0.77351965, 0.5470393 , 1. ]],
[[ 1. , 1. , 1. ],
[ 1. , 1. , 1. ],
[ 1. , 1. , 1. ]]])
Computes the r*e^(I*theta) form of derivatives from the grid of points. The derivatives are computed using quick-and-dirty taylor expansion and assuming analyticity. As such get_derivatives is primarily intended to be used for comparisions in plot_spiderweb and not for applications that require great precision.
INPUT:
OUTPUT:
EXAMPLES:
Standard usage with compute_on_grid:
sage: from sage.calculus.riemann import get_derivatives
sage: f(t) = e^(I*t) - 0.5*e^(-I*t)
sage: fprime(t) = I*e^(I*t) + 0.5*I*e^(-I*t)
sage: m = Riemann_Map([f], [fprime], 0)
sage: data = m.compute_on_grid([],19)
sage: xstep = (data[2]-data[1])/19
sage: ystep = (data[4]-data[3])/19
sage: dr, dtheta = get_derivatives(data[0],xstep,ystep)
sage: dr[8,8]
0.241...
sage: dtheta[5,5]
5.907...