3
$\begingroup$

Michael Trott in his book The Mathematica GuideBook for Symbolics, p. 1004, has illustrated a nice way to visualize the Riemann surface of the incomplete gamma function $\Gamma(\alpha, z)$. To plot the Riemann surface of, for example, $\Gamma\left(\frac{3+i}{10}, z\right)$, we have:

With[{α = 0.3 + 0.1 I, ε = 10^-12},
     GraphicsRow[Graphics3D[{EdgeForm[Thickness[0.002]],
                             SurfaceColor[Hue[0.09], Hue[0.18], 2.3],
                             Table[Last /@ Partition[Cases[
                 ParametricPlot3D[{r Cos[φ], r Sin[φ],
                                   #[Exp[2 k π I α] Gamma[α, r Exp[I φ]] +
                                   (1 - Exp[2 k π I α]) Gamma[α]]},
                                  {r, 0, 2}, {φ, -π + ε, π - ε},
                                  PlotPoints -> {30, 40}], _Polygon, ∞], 2],
                                   {k, -2, 2}]}, BoxRatios -> {1, 1, 2.5},
                            PlotRange -> All] & /@ {Re, Im}]]

where we have used the identity:

$$\Gamma\left(\alpha, \exp(2 k \pi i)z\right) = \exp(2 k \pi i \alpha)\, \Gamma\left(\alpha, z\right) + (1- \exp(2 k \pi i \alpha)) \Gamma\left(\alpha\right)$$

However, Mathematica does not return anything.

Any help is appreciated!

$\endgroup$
3
$\begingroup$

I believe that code is obsolete and was written before the use of GraphicsComplex. It could be updated, but perhaps the easiest way to solve your problem is to simply remove the Cases operator (and other constructs) which incorrectly extracts the polygon data and just plot the 5 sheets of the function using the code below. Note that I explicitly extract the real part of the function with the Re operator. Can plot the imaginary with Im. You can then experiment with the code. For example, plot one sheet, let $k=0$ for example, or plot the sheets with different colors, etc.

α = 0.3 + 0.1 I;
ee = 10^(-12)
myTable = 
  Table[
    ParametricPlot3D[
      {r Cos[φ], 
       r Sin[φ], 
       Re[Exp[2 k Pi I α] Gamma[α, r Exp[I φ]] + (1 - Exp[2 k Pi I α]) Gamma[α]]}, 
      {r, 0, 2}, {φ, -Pi + ee, Pi - ee}, 
      PlotPoints -> {30, 40}, 
      BoxRatios -> {1, 1, 1}], {k, -2, 2}]; 
Show[myTable]

enter image description here

$\endgroup$
3
  • $\begingroup$ Thanks for your answer. Instead of using the the mentioned identity, if we write the code as: $\endgroup$
    – user67559
    Sep 24 '19 at 22:29
  • $\begingroup$ g[z_, n1_] := Gamma[0.3 + 0.1 I, Abs[z] ((E^(I Arg[z])) E^(I 2 n1 Pi))]; p1 = ParametricPlot3D[{Re[z], Im[z], Im[g[z, 0]]} /. z -> r Exp[I t], {r, 0, 3}, {t, -Pi, Pi}, PlotStyle -> Red]; p2 = ParametricPlot3D[{Re[z], Im[z], Im[g[z, -1]]} /. z -> r Exp[I t], {r, 0, 3}, {t, -Pi, Pi}, PlotStyle -> Green]; Show[{p1, p2}, PlotRange -> All, BoxRatios -> {1, 1, 1}] $\endgroup$
    – user67559
    Sep 24 '19 at 22:29
  • $\begingroup$ It seems Mathematica produces the same plot for the first two sheets: n1=0 and n1=-1. Is my code correct? $\endgroup$
    – user67559
    Sep 24 '19 at 22:30
2
$\begingroup$

It seems to me that to fully reproduce what Michael Trott's code does, we have to write something like the following which shows both the real and the imaginary surfaces.

α = 0.3 + 0.1 I;
ee = 10^(-12);
indxs = Range[-2, 2];
colors = Most[Hue /@ Subdivide[Length[indxs]]];
plotF =
  Function[u,
    MapThread[
      ParametricPlot3D[
        {r Cos[φ],
         r Sin[φ],
         u[Exp[2 #1 Pi I α] Gamma[α, r Exp[I φ]] + (1 - Exp[2 #1 Pi I α]) Gamma[α]]}, 
        {r, 0, 2}, {φ, -Pi + ee, Pi - ee},
        PlotPoints -> {30, 40},
        BoxRatios -> {1, 1, 1},
        PlotStyle -> #2] &, 
      {indxs, colors}]];
Column[
  Legended[
    Show[plotF[#], PlotRange -> All, ImageSize -> 400],
    SwatchLegend[colors, indxs, LegendLabel -> Style["k", "TBI", 16]]] & /@ {Re, Im}]

plots

$\endgroup$
1
  • $\begingroup$ Many thanks for your answer. $\endgroup$
    – user67559
    Sep 25 '19 at 9:14
2
$\begingroup$

For reference, this is what the code produces in version 5.2 (the latest version at the time the bulk of Trott's book was written):

version 5.2 screenshot

Now, the reason for the Cases[(* stuff *), _Polygon, ∞] and the Last /@ Partition[(* stuff *)] business is that it was Trott's method to help reveal the structure of the Riemann surface, by extracting all polygons from the plots, and then deleting every other polygon in the plotted surfaces.

In the current version of Mathematica, a little more work is necessary to reproduce this figure. Instead of manually deleting polygons, we can instead fiddle with the Mesh and MeshShading settings, like so:

With[{α = 0.3 + 0.1 I, ε = 1*^-12},
     GraphicsRow[Table[
                 ParametricPlot3D[Table[{r Cos[φ], r Sin[φ],
                                         ff[Exp[2 k π I α] Gamma[α, r Exp[I φ]] +
                                            (1 - Exp[2 k π I α]) Gamma[α]]},
                                        {k, -2, 2}],
                                  {r, 0, 2}, {φ, -π + ε, π - ε},
                                  BoxRatios -> {1, 1, 2.5}, Lighting -> "Classic", 
                                  Mesh -> {31, 41},
                                  MeshShading -> {{None, Automatic}, {Automatic, None}}, 
                                  MeshStyle -> Thickness[0.002], PlotRange -> All,
                                  PlotStyle -> Directive[Hue[0.09],
                                                         Specularity[Hue[0.18], 2.3]]],
                            {ff, {Re, Im}}]]]

Riemann surface of upper incomplete gamma function, cutaway


Nowadays, however, we have at our disposal a graphics directive that Trott didn't: Opacity[]. We can use it to have a good look at the Riemann surface's structure by making the surface translucent, and forgo all the polygon removal business:

With[{α = 0.3 + 0.1 I, ε = 1*^-12},
     GraphicsRow[Table[ParametricPlot3D[Table[{r Cos[φ], r Sin[φ],
                                               ff[Exp[2 k π I α] Gamma[α, r Exp[I φ]] +
                                               (1 - Exp[2 k π I α]) Gamma[α]]},
                                              {k, -2, 2}],
                                        {r, 0, 2}, {φ, -π + ε, π - ε},
                                        BoxRatios -> {1, 1, 2.5}, Mesh -> None,
                                        PlotRange -> All,
                                        PlotStyle -> Directive[Hue[0.09], Opacity[0.6]]],
                       {ff, {Re, Im}}]]]

Riemann surface of upper incomplete gamma function, translucent

$\endgroup$
0
$\begingroup$

Here is code to plot two sheet separately. You can add to it:

k = 0; α = 0.3 + 0.1 I; ε = 10^(-12)
p1 = ParametricPlot3D[{r Cos[φ], r Sin[φ], 
                       Re[Exp[2 k π I α] Gamma[α, r Exp[I φ]] +
                          (1 - Exp[2 k π I α]) Gamma[α]]},
                      {r, 0, 2}, {φ, -π + ε, π - ε}, PlotPoints -> {30, 40}, 
                      BoxRatios -> {1, 1, 1}, PlotStyle -> Red];
k = -1;
p2 = ParametricPlot3D[{r Cos[φ], r Sin[φ], 
                       Re[Exp[2 k π I α] Gamma[α, r Exp[I φ]] +
                          (1 - Exp[2 k π I α]) Gamma[α]]},
                      {r, 0, 2}, {φ, -π + ε, π - ε}, PlotPoints -> {30, 40}, 
                      BoxRatios -> {1, 1, 1}, PlotStyle -> Blue];
Show[{p1, p2}, PlotRange -> All]

Riemann surface

$\endgroup$
1
  • $\begingroup$ Another question. By analytic continuation we can ascertain that the identity (regarding the gamma function) used in the code also holds when alpha is zero. However, if we set alpha = 0, we get an error. How can we fix that? $\endgroup$
    – user67559
    Sep 25 '19 at 2:06

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy