Plotting graphs of Modular Forms

Question

After watching all the 8 parts of "“Introduction to Modular Forms,” by Keith Conrad" on YouTube, I got "extremely intrigued" by plotting graphs of Modular Forms ( on SL(2,Z) ).

So after watching all those videos, I tried the following approach by dipping my hands into SageMath fundamentally by obtaining the q-expansion and using the Fourier series to reconstruct it.

I used the following "procedure" ( this has been done using SageMath ) :

1. obtain the modular form M4 = ModularForms(SL2Z,4)

2. obtain 40 coefficients of the q-expansion print (M4.q_expansion_basis(40))

3. define FM(q) = 1 + 240q + 2160q^2 + 6720q^3 + 17520q^4 + 30240q^5 + 60480q^6 + 82560q^7 + 140400q^8 + 181680q^9 + 272160q^10 + 319680q^11 + 490560q^12 + 527520q^13 + ..etc. etc. .. up to .. +14770560q^39

4. define FQ(z) = e^(2Ipi*z)

5. finally plot the modulus of FM(FQ(z)) with a colormap ( i.e. complex_plot(FM4(FQ()), (-1,1),(0,1.5)))

I get results like this picture :

But as you can see "the image definition is not the best" , I was looking around for papers/things about "plotting modular forms" could not find much, except one paper mentioning something like "using 5000 terms of the Fourier Serie".

Now the question is "in any better way to get more detailed pictures of Modular forms ?", I suspect mine is a bit of a "brute force approach" also probably prone to much numerical instability ?

If you ask me "why all that ?" .. "because I wanted to see how they look like" :)

Many thanks in advance for any reply.

Cheers.

Answer 1

You might be interested in David Lowry-Duda's work on that.

• David Lowry-Duda. Visualizing modular forms. arXiv:2002.05234

• David Lowry-Duda. Notes behind a talk: visualizing modular forms. Blog post, 2019-11-22. https://davidlowryduda.com/notes-behind-a-talk-visualizing-modular-forms/

• David Lowry-Duda. phase_mag_plot, a Sage package for plotting complex functions.

  • blog post: https://davidlowryduda.com/phase_mag_plot-a-sage-package-for-plotting-complex-functions/
  • code repository: https://github.com/davidlowryduda/phase_mag_plot/

Answer 2

There are two things necessary to increase the definition of the image:

1. Compute the modular form to higher accuracy, and
2. Plot the modular form in a finer mesh.

In sage, the way to do the latter is with the plot_points parameter in complex_plot. By default, sage will evaluate the function on a $100 \times 100$ grid and interpolate between the points. This gives a pretty low-definition picture.

Let's include runnable code. Although I separate the code into chunks, think of all the code as a single code file, in that later code snippets depend on earlier code snippets.

The following code redoes what you've done. It computes a degree $39$ polynomial using the first 40 coefficients of $f$, the unique modular form (actually an Eisenstein series) of weight $4$ on $\mathrm{SL}(2, \mathbb{Z})$. Then it plots it.

# sage code

Htoq = lambda x: exp(2 * CDF.pi() * CDF.0 * x)
DtoH = lambda x: (-CDF.0 * x + 1) / (x - CDF.0)
Dtoq = lambda x: Htoq(DtoH(CDF(x)))
C.<t> = CC[]

M4 = ModularForms(1, 4)
f = M4.basis()[0]
coeffs = f.coefficients(list(range(40)))  # take 40 coefficients
fpoly = C(coeffs)  # interpret coeffs in the polynomial ring C.<t>

P = complex_plot(
  lambda z: fpoly(Htoq(z)), 
  (-1, 1), (0.01, 2),
  aspect_ratio=1, figsize=[4, 4])
P.axes(show=False)
P

This makes the following image, from $40$ coefficients evaluated on a $100 \times 100$ grid.

We can instead evaluate it on a $500 \times 500$ grid with the following code change.

# continuing from above code
P = complex_plot(
  lambda z: fpoly(Htoq(z)), 
  (-1, 1), (0.01, 2),
  aspect_ratio=1, figsize=[4, 4],
  plot_points=500   # <--- new line
)
P.axes(show=False)
P

Notice the higher resolution. The bands of color at the bottom are numerical artifacts from low precision evaluation of the modular form.

Using more coefficients is one way to increase the accuracy near the boundary. See the following code.

coeffs = f.coefficients(list(range(100)))  # <-- 100 coefficients
fpoly = C(coeffs)

P = complex_plot(
  lambda z: fpoly(Htoq(z)), 
  (-1, 1), (0.01, 2),
  aspect_ratio=1, figsize=[4, 4],
  plot_points=500)
P.axes(show=False)
P

This gives the following picture.

Now there is almost no banding at the boundary. Similar techniques would work for other visualizations, such as on the Poincare disk.

P = complex_plot(
  lambda z: +Infinity if abs(z) >= 0.99 
            else fpoly(Dtoq(z)) * exp(1.2 * CDF.pi() * CDF.0),
  (-1, 1), (-1, 1), aspect_ratio=1, figsize=[4,4],
  plot_points=500)
P.axes(show=False)
P

This gives the following image on the disk.

There is some banding around parts of the boundary from approximation problems.

I note that the argument factor exp(1.2 * pi * i) is there just to make the phase corresponding to real values blue, which is more pleasant than the default red. For the older images corresponding to modular forms in the LMFDB, we did the same thing more or less. (We now produce images like on this page, which use a different colormap. Reasons and methods for this are noted in this paper).

I'll also note that there are other ways to compute the modular form to higher accuracy. For level $1$, you can pullback each point to the fundamental domain (where the imaginary part is large, and thus the $q$-expansion converges rapidly). It is necessary to track the matrices used and use the modularity $f(\gamma z) = (cz + d)^k f(z)$.

I'll give sample code that does this, but with a small caveat: I use this code outside of sage in ways that allow this to be done faster or in different visualizations. I put code modified from what was used to make this video. It would be possible to do this a bit more simply using sage's concept of matrices, matrix multplication, and actions.

# sage code

Htoq = lambda x: exp(2 * CDF.pi() * CDF.0 * x)
DtoH = lambda x: (-CDF.0 * x + 1) / (x - CDF.0)
Dtoq = lambda x: Htoq(DtoH(CDF(x)))
C.<t> = CC[]

M4 = ModularForms(1, 4)
f = M4.basis()[0]
coeffs = f.coefficients(list(range(40)))  # only 40 coefficients
fpoly = C(coeffs)

def in_fund_domain(z):
    x = z.real()
    y = z.imag()
    if x < -0.51 or x > 0.51:
        return False
    if x*x + y*y < 0.99:
        return False
    return True

def act(gamma, z):
    a, b, c, d = gamma
    return (a*z + b) / (c*z + d)

def mult_matrices(mat1, mat2):
    a, b, c, d = mat1
    A, B, C, D = mat2
    return [a*A + b*C, a*B + b*D, c*A + d*C, c*B + d*D]

Id = [1, 0, 0, 1]

def pullback(z):
    """
    Returns gamma, w such that gamma(z) = w and w is
    (essentially) in the fundamental domain.
    """
    z = CDF(z)
    gamma = Id
    count = 1
    while not in_fund_domain(z):
        if (count > 1000): print(z, count)
        count += 1
        x, y = z.real(), z.imag()
        xshift = -floor(x + 0.5)
        shiftmatrix = [1, xshift, 0, 1]
        gamma = mult_matrices(shiftmatrix, gamma)
        z = act(shiftmatrix, z)
        if x*x + y*y < 0.99:
            z = -1/z
            gamma = mult_matrices([0, -1, 1, 0], gamma)
    return gamma, z

def smart_compute(z):
    gamma, z = pullback(DtoH(z))
    a, b, c, d = gamma
    return (-c*z + a)**4 * fpoly(Htoq(z))

P = complex_plot(
  lambda z: +Infinity if abs(z) >= 0.99
            else smart_compute(z) * exp(1.2 * CDF.pi() * CDF.0),
  (-1, 1), (-1, 1), aspect_ratio=1, figsize=[6, 6],
  plot_points=500)

P.axes(show=False)
P

Using only 40 coefficients, this gets a very high fidelity image (but takes approximately 20 seconds to run on my machine).

For higher level, it is not possible to always choose points in the fundamental domain with large imaginary part. Alternative methods are available, but they're frequently rather annoying in my experience.