I'm considering solving elliptic functions over a fundamental domain of the torus with half-periods $\omega_{1}=\pi/2$ and $\omega_{2} = \pi \tau /2$, where $\tau$ is the modular parameter of the torus. The equation I want solutions of is: $$\wp(u \, | \, \omega_{1}, \omega_{2})= -\frac{1}{3}E_{2}(\tau)$$

which definitely must have two solutions $u$ in the parallelogram for all $\tau$. Or perhaps one solution or order 2.

My main issues are not knowing how to effectively parameterize the fundamental parallelogram in Mathematica, as well as my NSolve routine not working properly. My idea was to fix $\omega_{1}$ and $\omega_{2}$ and then have mathematica consider the domain: $$\{x\omega_{2} + y \omega_{2} \, | \, 0 \leq x \leq 2, -1 \leq y \leq 1\} \subseteq \mathbb{C}.$$

I think this is OK, but I'm also worried it might be a bad parameterization causing my NSolve issues. The code I have is:

tau = 0 + (3/2)*I; w1 = Pi/2; w2 = Pi*tau/2; inv = WeierstrassInvariants[{w1, w2}]; E2[t_] = 1 - 24*Sum[(n*Exp[2*Pi*I*(t)*n])/(1 - Exp[2*Pi*I*(t)*n]), {n, 1, 300}]; WP[x_, y_] = WeierstrassP[w1*x + w2*y, inv]; L = -(1/3)*N[E2[tau], 50]; NSolve[{WP[x, y] == L && 0 <= x <= 2 && -1 <= y <= 1}, {x, y}, WorkingPrecision -> 50]

The problem is that NSolve is just spitting my expression back out immediately without computing anything. I've tried a number of things including telling Mathematica to do it over the Reals, using a single complex variable $u$ instead of the $x,y$, as well as dropping the domain restrictions, and nothing seems to work.

My function E2[t_] for the Eisenstein series isn't the problem; it spits out a number incredibly quickly.

I'd really like to avoid using FindRoot if at all possible because I don't want to have to estimate where these two special points are in the parallelogram.

Thanks a lot for any advice!

  • 1
    $\begingroup$ Have you tried making a plot? $\endgroup$ Dec 11, 2015 at 20:59
  • $\begingroup$ @J.M. So I tried using the function about halfway down that link by Vitaliy Kaurov. That example worked perfectly, but when I applied it to my case I got a long list of wild errors that I've never seen in Mathematica. I feel like that code was pretty transparent but something didn't work. I used that code twice for the real and imaginary parts of my WP function above. Is there another of the codes in that link that you like? $\endgroup$
    – Benighted
    Dec 11, 2015 at 21:25
  • 1
    $\begingroup$ What I had in mind was for you to consider Re[WP[x, y] - L] == 0 and Im[WP[x, y] - L] == 0 as simultaneous equations that you could feed into any of the answers to the linked question. $\endgroup$ Dec 11, 2015 at 21:27
  • $\begingroup$ That worked great, thanks! Produced unexpected results, but I think the plots themselves are working. $\endgroup$
    – Benighted
    Dec 11, 2015 at 22:23

1 Answer 1


Here's a polynomial interpolation method, which can be be found in Chapter 5 of Boyd (2014).

nn = 64;
z0 = w1 + w2;
rr = 1.1 w1;
ff = N[WeierstrassP[z0 + rr #, inv] - L, Precision[#]] &;
wprec = MachinePrecision;
tj = 2 Pi*Range[0, nn - 1]/nn;
wj = N[Exp[I tj], wprec];
fj = ff /@ wj; (* f[zj] *)
aa = InverseFourier[fj]/Sqrt[nn];

(* Rough check of accuracy of interpolation *)
"condition"@# -> Log10@Ratios[#] &@ N@ MinMax@ Abs@ fj
ip = FromDigits[Reverse@aa, (z - z0)/rr];
Max@ Table[(WeierstrassP[z, inv] - L - ip)/(0WeierstrassP[z, inv] - L) /. 
    z -> z0 + r Exp[I t] // N // Abs,
  {r, (1/8 - 0.01) rr, rr, rr/8}, {t, 0., 2 Pi, 0.2/r}]
  "condition"[{0.0694093, 0.75655}] -> {1.03742}

z2 = Eigenvalues@ companionMatrix[aa];
roots = z0 + rr*Select[z2, Abs[#] < 0.999 &]
(*  {0.776416 + 2.35619 I, 2.36518 + 2.35619 I}  *)

  EdgeForm@Gray, LightBlue, Disk[ReIm@z0, rr], Gray, Point[ReIm@z0],
  Black, PointSize[Medium], Point@ReIm[roots]},
 Frame -> True]

Mathematica graphics

Auxiliary code:

For the companion matrix, you can use

companionMatrix[coeffs_] := 
    Band[{1, 2}] -> 1, {Length@coeffs - 2, Length@coeffs - 1}],
   {-coeffs[[;; -2]]/coeffs[[-1]]}


companionMatrix = NRoots`CompanionMatrix
  • $\begingroup$ You've taken the polynomial route way further than I have. Awesome. :) $\endgroup$ May 5, 2016 at 14:46
  • $\begingroup$ @J.M. Thanks! I may add some explanation when I'm done grading final exams. (You have to avoid the poles, for instance.) I taught some of the real function, colleague matrix stuff this semester (my first time teaching intro. num. anal.). I have you to thank for the Boyd reference. :) $\endgroup$
    – Michael E2
    May 5, 2016 at 15:00

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.