Let us demonstrate the steps outlined above using the following not-too-simple function.

A contour plot shows that wExample has a complicated branch cut structure. As mentioned in the previous installment, in a contour plot the branch cuts show up as clusters of contour lines.

We canonicalize all functions by writing them explicitly in logarithms. This rewriting allows for a more uniform treatment of all functions. We will not have to deal separately with all of the individual Arc* functions. At the same time, we do not want meromorphic functions like Sin to be rewritten as a sum of exponentials, because potentially this could multiply the number of logarithmic terms. Instead we apply the function TrigToExp only to the Arc* functions.

Here are two examples of canonicalized functions.

This representation shows the source of the three poles visible in the above contour plot. They are caused by the arguments of the logarithms when the arguments become zero. This event happens at exactly three -values.

The next step is to determine all potential branch points. (I say potential branch points because in some situations they might turn out not to be branch points, but including them does not hurt in any way.) To do this we select all -dependent log terms and power terms (with fractional exponents). The function multiValuedTerms returns a list of all multivalued parts of the function under consideration. We are looking for the multivalued functions of the form , , , and Root[,].

Here are two examples of the result of applying multiValuedTerms.

Here are the multivalued terms of our example function.

The branch points are then given as the zeros and infinities of these expressions, because a rational function can map the branch point at infinity to a finite -value. (An example of this situation is , where at the -branch point of the square root function reappears.) We equate the denominator and the numerator independently to 0. When calling Solve, error messages may be generated (either because the equations could not be solved at all or because not all branch points could be found). If messages were generated, we give up.

The branch points of are the values of such that or .

The branch points of and are the values of such that or .

The branch points of Root[, ] are the values of such that has a multiple root with respect to . (Here is a pure function.)

Here are some examples of branchPoints at work. In the first example the branch point originates from the Log term, and the other five branch points come from the Root-object.

We encourage the reader to carefully examine the following values for the potential branch points and to note from which functions they arise.

The following function is meromorphic in and does not have any branch points.

We get the following branch points for our example function wExample.

We transform all exact numeric quantities into approximate floating point numbers. All of the further work is mainly numerical in nature anyway.

Now that we have found the finite branch points, we divide the complex -plane into sectors between them by using a tensor product structure formed by the radial and azimuthal positions of the branch points. We use a polar coordinate system for doing this. (Of course, one could also use a Cartesian one--we leave this to the interested reader.) The sectors returned from sectors are strictly inside the areas bounded by either the potential branch points (achieved by adding/subtracting a small to all coordinate values at branch points) or regular points. For a graphical representation of the Riemann surface, which is the only goal here, one will usually not see the resulting little gaps.

These are the sectors originating from the example function.

Let us visualize the sectors together with the branch points by coloring each sector with a randomly chosen color. The function sectorPolygon generates 2D polygons representing the individual sectors.

Now we must consider all the multivalued terms in our function to properly continue them analytically. Be aware that the sector boundaries have nothing in common with the branch cuts of the original functions; only the branch points coincide. The actual continuation will be carried out in this way.

a. Radicals will be continued according to , .

b. Logarithms will be continued according to , . We will use only a finite number of sheets--typically three--later on for the graphics.

c. Powers of the form will be continued by rewriting them using logarithms as . This reduces this case to case b.

d. Root objects will be continued according to Root[,]Root[,], .

To differentiate the various multivalued terms, we give each of them an "index," like or . This index allows us to keep track of them uniquely later on. This is important because equal-looking terms could give rise to different sheets, like in the example . The two square roots are unrelated and each one generates two branches, so that the resulting Riemann surface has four sheets.

Here are all multivalued terms numbered in wExample.

The numbering for two other examples is as follows.

Now we extract all indexed terms from the function and collect them in a list.

Each of the numbered multivalued functions is now replaced by a symbolic function . After rewriting them into a nicer form, we can deal with them numerically later on inside NDSolve.

The result of substitutions is a list with three elements. The first element is a list of the replacements of the multivalued functions by the . The second element is a similar list without the recursive back-substitution of the . The third element is the original function expressed through the s. For our example function application of substitutions results in the following six substitutions.

We now generate differential equations from the last list of substitutions. Doing this seems like a lot of work, considering that we already know the solution of these differential equations--these are just the values at the other sheets--but it is unavoidable. Although we know the location of the branch points, we do not know the detailed location of the branch cuts. So our solution might not represent the intended sheet. The branch cuts might form complicated patterns inside the sectors calculated above, and it would be very difficult or virtually impossible to construct the "correct" sheet purely symbolically. (Consider log(aComplicatedPolynomial). The branch points are just the zeros of aComplicatedPolynomial, but the branch cuts will typically wind through the individual sectors as complicated curves.) The resulting system of differential equations is then of the form where all of the are now meromorphic in all arguments, and no multivaluedness is present anymore in the solution of the differential equation system for the given initial conditions.

We obtain the differential equation for by using straightforward differentiation.

We use implicit differentiation to obtain the differential equation, meromorphic in and for .

We obtain the differential equation for Root[,] by implicit differentiation of . Again is a pure function. Observe the parentheses in . First we differentiate the pure function ξ which parametrically depends on z with respect to z, and then we apply the resulting pure function to the argument w.

Here are the differential equations corresponding to subs.

Here the differential equation corresponding to a Root-object is calculated.

The right-hand sides of the last system of differential equations often still contain the derivatives of some of the . We recursively insert the corresponding right-hand sides to get a system of differential equations with rational right-hand sides that can be solved directly by NDSolve.

Let us look at the resulting system of differential equations for two more examples. Here is the differential equation system corresponding to . The function generates the system of differential equations in the form used to visualize Riemann surfaces.

Here are the systems of differential equations that corresponds to some more complicated functions.

Now we must consider the sheets we actually want to display. When near a branch point (i.e., at the corners of our sectors), we will choose the sheets that we want to continue analytically. For radicals we will use all possible sheets (if the reader is going to use high order radicals one might want to restrict this). For logarithms we will use by default three neighboring sheets (the one "below" and the one "above" the "original" one). As with radicals, we take all sheets into account when using Root-objects.

This method results in the following sheet realizations of our numbered functions.

Forming the outer product of all possible sheet realizations of all multivalued functions, we get all potential sheets.

For our example function there are potential sheets. Usually, however, we will not really have so many different sheets! When we calculate numerical values of the sheets (see below), some of the sheets might turn out to be numerically the same. A simple example is This expression contains one square root to be continued, but the canonicalized form contains two copies of the same square root which are differently numbered. Later on we will discuss how to avoid displaying the same sheet twice.

Here is the first realization.

Now we implement two functions that transform a given differential equation of the form into one in which or ( or being polar coordinates) stay constant. The dependent variable remains unchanged.

We get the following differential equations for our example function.

Now we will finally solve the coupled set of differential equations to get the Riemann surface. The function sheetSectorPatch generates the polygons for one sheet realization within one sector.

We set the MaxSteps, PrecisionGoal, and AccuracyGoal options of NDSolve to values higher than the default. This makes sure that the examples treated in the next section will work fine. For more complicated functions, it may be necessary to use more steps and a higher precision. The use of high-precision arithmetic via WorkingPrecision->valueLargerThan$MachinePrecision might also be needed.

The following shows one sheet inside the first sector for wExample.

One sheet (not always the same) for all sectors is shown.

Here is a view of the surface patches from the last picture from "above." The sector structure is clearly visible.

We now consider the generation of all sheets of interest inside one sector. The first step is to calculate the numerical values of the sheets at one of the corners of a sector. If two such values agree, then they represent the same sheet (both are generated by a differential equation without multiple solutions for given initial conditions). The built-in function Union relies on sorting and is not applicable for unioning complex numbers that do not form a well-ordered set. The following function complexUnion implements the needed unioning for the complex starting values. Be aware that the function complexUnion has complexity in comparison to the much smaller complexity of Union--but for the short lists we have here, this speed/complexity difference does not matter.

The function differentSheetRealizations returns a list of all the different initial conditions present in at the point .

Now we generate all patches inside a given sector. After we apply Re or Im or the like to the solutions of the differential equations, some sheets might coincide. It is not a good thing to have a graphics expression containing nearly identical polygons because they intersect frequently and a lot of time is spent in calculating and rendering their intersections. To avoid this situation we implement a function dropMultipleSheets which eliminates such multiple polygons/sheets.

The function sectorSheetPatches finally produces all sheets as specified by sheetRealizations.

Here all the patches for our example function for the fifth sector are shown. One easily sees the branch point where the sheets touch each other.

The function patchPlotPoints divides the plot points in radial and azimuthal directions in a "fair" manner between the individual sectors.

Now we finally have all the ingredients together to implement the function RiemannSurface which will generate the polygons of a few sheets of the Riemann surface of the function w. We add a few simple tests to keep inappropriate functions from getting passed along.

The function notTreatedFunctions returns a list of all functions that are not dealt with in this implementation of construction of Riemann surfaces. foreignFunction collects all functions that are not explicitly listed as "known" functions.

The basic test numericQ checks to see if w is a numeric function of the independent variable z at all points (to catch inputs similar to RiemannSurface[{Sqrt[x]+Log[x],z},...]).

The function RiemannSurface integrates all of the above step-by-step generation of the sheets into one function.

For convenience, we now allow equations as inputs so that implicit forms can be used, too.

Now we have implemented all functions needed to visualize the Riemann surface. For our example function we arrive at the following pictures for the imaginary and real parts.

As the previous picture shows, the imaginary part is a faithful image of the Riemann surface.

Here is another example: the Riemann surface of an algebraic function.

Note: we can get a purely polynomial description of this function by using the function Resultant.

Here is another example that contains the inverse of the error function.

You did not really expect the last example to work, did you? As mentioned in the introduction, the Riemann surface of such functions will be constructed in Part III of this series, and this is only the third installment of Part II! According to the test implemented in notTreatedFunctions the foreignFunction InverseErf is caught and detected.