Here is a computational problem given to me by Vinayak Vatsal. Given the Riemann theta function

and related Riemann theta function with characteristic $[\alpha, \beta]$

compute $\theta[\sigma,0](0, 2\tau)$ for $\tau = zA, zB$ and various $z \in \mathbb{C}, \text{Im}(z) > 0$ with $\sigma \in {0, 1/2}^3$. The definitions of the matrices $A,B \in \mathbb{C}^{3 \times 3}$ are

Below are some computations done in an iPython notebook using abelfunctions. You can download the notebook itself here, which requires that abelfunctions be installed in order to execute.

The Setup

In [1]:

import numpy
import matplotlib
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1 import ImageGrid
import matplotlib.colors

from abelfunctions import RiemannTheta

In [2]:

A = numpy.matrix(
    [[2, 1, 1],
     [1, 4, 1],
     [1, 1, 6]],
    dtype = numpy.complex)
B = numpy.matrix(
    [[23, -5, -3],
     [-5, 11, -1],
     [-3, -1, 7]],
    dtype = numpy.complex)

In [3]:

# construct the funciton we want to compute the values of:
#
# theta[sigma,0](0,zM)
#
# where sigma, z, and M can vary

def vitsal(z, sigma, mat):
    # construct the Riemann (period) matrix tau and 2tau
    tau = z*mat
    tau2 = 2*tau

    # compute theta with characteristic (RiemannTheta.characteristic()
    # is untested.)
    dot = numpy.dot(sigma.T, numpy.dot(tau2, sigma))
    scale = numpy.exp(2*numpy.pi*1.0j * dot).item()
    arg = numpy.dot(tau2, sigma)
    return scale * RiemannTheta(arg, tau2)

# vectorize over z since we want to see the behavior of
# this function for some varying z
vitsal = numpy.vectorize(vitsal, excluded=['sigma', 'mat'])

Fixed $z \in \mathbb{C}$, All Characteristics $\sigma \in {0,1/2}$

We (arbitrarily) choose $z = i$ and $z = 2i$. Although, the plots in the following section suggests that the choice of $z$ in the upper half plane doesn’t matter except, possibly, in convergence rate of the implementation of the Riemann theta function to the suggested value of the function. That is, it seems like the truncated sum of the Riemann theta function converges to the actual value as we add more terms faster for $\text{Im}(z)$ large.

In [4]:

def print_values(z, mat):
    # constuct list of all half integer characteristics
    halfs = [0,0.5]
    sigmas = []
    for s1 in halfs:
        for s2 in halfs:
            for s3 in halfs:
                sigmas.append(numpy.matrix([s1,s2,s3]).T)

    # for each characteristic, print characteristic and
    # value of theta function
    for sigma in sigmas:
        val = vitsal(z, sigma=sigma, mat=mat)
        print 'sigma.T: %s\tvalue: %s'%(sigma.T, val)

In [5]:

print_values(1.0j, A)

Out[5]:

sigma.T: [[0 0 0]]      value: (1.00000697473+0j)
sigma.T: [[ 0.   0.   0.5]]      value: (6.52457377777e-09+0j)
sigma.T: [[ 0.   0.5  0. ]]      value: (6.9877095594e-06+0j)
sigma.T: [[ 0.   0.5  0.5]]      value: (2.2795992109e-14+0j)
sigma.T: [[ 0.5  0.   0. ]]      value: (0.00373488548776+0j)
sigma.T: [[ 0.5  0.   0.5]]      value: (1.21842677994e-11+0j)
sigma.T: [[ 0.5  0.5  0. ]]      value: (1.30369858713e-08+0j)
sigma.T: [[ 0.5  0.5  0.5]]      value: (1.58402415218e-19+0j)

In [6]:

print_values(2.0j, A)

Out[6]:

sigma.T: [[0 0 0]]      value: (1.00000000002+0j)
sigma.T: [[ 0.   0.   0.5]]      value: (4.24116597336e-17+0j)
sigma.T: [[ 0.   0.5  0. ]]      value: (1.21615991209e-11+0j)
sigma.T: [[ 0.   0.5  0.5]]      value: (5.15790006267e-28+0j)
sigma.T: [[ 0.5  0.   0. ]]      value: (6.97468471242e-06+0j)
sigma.T: [[ 0.5  0.   0.5]]      value: (1.47903977386e-22+0j)
sigma.T: [[ 0.5  0.5  0. ]]      value: (4.24116597336e-17+0j)
sigma.T: [[ 0.5  0.5  0.5]]      value: (6.27280941128e-39+0j)

In [7]:

print_values(1.0j, B)

Out[7]:

sigma.T: [[0 0 0]]      value: (1+0j)
sigma.T: [[ 0.   0.   0.5]]      value: (2.81426845749e-10+0j)
sigma.T: [[ 0.   0.5  0. ]]      value: (9.81431759353e-16+0j)
sigma.T: [[ 0.   0.5  0.5]]      value: (1.47903461596e-22+0j)
sigma.T: [[ 0.5  0.   0. ]]      value: (4.16240046723e-32+0j)
sigma.T: [[ 0.5  0.   0.5]]      value: (1.79873633572e-33+0j)
sigma.T: [[ 0.5  0.5  0. ]]      value: (1.79873633572e-33+0j)
sigma.T: [[ 0.5  0.5  0.5]]      value: (4.16240046723e-32+0j)

In [8]:

print_values(2.0j, B)

Out[8]:

sigma.T: [[0 0 0]]      value: (1+0j)
sigma.T: [[ 0.   0.   0.5]]      value: (7.9201069508e-20+0j)
sigma.T: [[ 0.   0.5  0. ]]      value: (9.63208298267e-31+0j)
sigma.T: [[ 0.   0.5  0.5]]      value: (2.18754339521e-44+0j)
sigma.T: [[ 0.5  0.   0. ]]      value: (1.73255776496e-63+0j)
sigma.T: [[ 0.5  0.   0.5]]      value: (3.23545240544e-66+0j)
sigma.T: [[ 0.5  0.5  0. ]]      value: (3.23545240544e-66+0j)
sigma.T: [[ 0.5  0.5  0.5]]      value: (1.73255776496e-63+0j)

Varying $z \in \mathfrak{h}_1$, $\sigma = [0,0,0]$ and $\sigma = [1/2,0,0]$

We now draw contour plots for the real, imaginary, and absolute value of the theta function for various $z$ in the upper half plane. We chose $\sigma = [0,0,0]$ because it was the only characteristic resulting in a theta value of 1. We chose $\sigma = [1/2,0,0]$ because it seems to be the one most sensitive to changing $z$ and the furthest away from zero for small values of $z$.

In [9]:

# computing the value of the theta function for various z

# set the sigma, z-grid size, and plotting window
sigma = numpy.matrix([0,0,0]).T
N = 64
zre_left, zre_right = -0.2,0.2
zim_left, zim_right = 1,1.2

# compute theta for all of these z
x = numpy.linspace(zre_left, zre_right, N)
y = numpy.linspace(zim_left, zim_right, N)
X,Y = numpy.meshgrid(x,y)
Z = X + 1.0j*Y
V = vitsal(Z, sigma=sigma, mat=A)

In [10]:

# a fancy plotting function so we can see the real, imaginary, and
# absolute values with the same colorbar. took more time to write this
# than the rest of the code in this notebook

def fancy_plot(X, Y, V, figwidth=10):
    # (... omitted for brevity ...)

In [11]:

fancy_plot(X, Y, V)

Out[11]:

In [12]:

# same as above, but with sigma = [0.5,0,0]

# set the sigma, z-grid size, and plotting window
sigma = numpy.matrix([0.5,0,0]).T
N = 64
zre_left, zre_right = -1,1
zim_left, zim_right = 1,2

# compute theta for all of these z
x = numpy.linspace(zre_left, zre_right, N)
y = numpy.linspace(zim_left, zim_right, N)
X,Y = numpy.meshgrid(x,y)
Z = X + 1.0j*Y
V = vitsal(Z, sigma=sigma, mat=A)

# plot
fancy_plot(X, Y, V)

Out[12]: