Using NSolve for Elliptic Equations over Fundamental Parallelogram in Complex Plane

Question

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!

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}
  1.7511*10^-12
*)

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

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

Auxiliary code:

For the companion matrix, you can use

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

Or

companionMatrix = NRoots`CompanionMatrix