# Using NSolve for Elliptic Equations over Fundamental Parallelogram in Complex Plane

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!

• Have you tried making a plot? – J. M. will be back soon Dec 11 '15 at 20:59
• @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? – Benighted Dec 11 '15 at 21:25
• 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. – J. M. will be back soon Dec 11 '15 at 21:27
• That worked great, thanks! Produced unexpected results, but I think the plots themselves are working. – Benighted Dec 11 '15 at 22:23

## 1 Answer

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 = NRootsCompanionMatrix
`
• You've taken the polynomial route way further than I have. Awesome. :) – J. M. will be back soon May 5 '16 at 14:46
• @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. :) – Michael E2 May 5 '16 at 15:00