Weierstrass P-function

[stla]

Hello,

It would be nice to have the Weierstrass P-function in mpmath.

[Hyrumveneno]

Agreed. Do you have a algorithm in mind to use?

[stla]

It is implemented in the R package elliptic. Does it help?

[fredrik-johansson]

It can be computed using Jacobi theta functions: https://fungrim.org/entry/af0dfc/

[Hyrumveneno]

It is implemented in the R package elliptic. Does it help?

Naw, I was just asking because if you know the algorithm, most of the time, it's not too hard to implement.

It can be computed using Jacobi theta functions: https://fungrim.org/entry/af0dfc/

It looks like you could implement it if you wanted to.

Sadly, I can't right now because I'm working on SymPy.

[stla]

I seriously believe this formula is wrong. This is much more complicated. This took me a long time but I've finally done an implementation. I compared the printed results with Wolfram and they are the same.

from mpmath import jtheta, pi, exp, sqrt, polyroots, agm

def p_weierstrass_from_g2_g3(g2, g3):
    r1, r2, r3 = polyroots([4, 0, -g2, -g3])
    e3 = r3
    a = sqrt(r1 - r3)
    b = sqrt(r1 - r2)
    c = sqrt(r2 - r3)
    if abs(a + b) < abs(a - b):
        b *= 1j
    if abs(a + c) < abs(a - c):
        c *= 1j
    if abs(c + 1j*b) < abs(c - 1j*b):
        e3 = r1
        a *= 1j 
        tmp = b
        b = 1j * c
        c = 1j * tmp
    w1 = 1 / agm(a, b)
    w3 = 1j / agm(a, c)
    q = exp(1j * pi * w3/w1)
    return lambda z: (
        e3 + (pi * jtheta(2,0,q) * jtheta(3,0,q) * jtheta(4,z/w1,q)
              / (pi * w1 * jtheta(1,z/w1,q)))**2    
    )

z = 1+1j
g2 = 5 - 3j
g3 = 2 - 7j
print(p_weierstrass_from_g2_g3(g2, g3)(z))
# (-0.0300998915997091 - 0.990427096602199j)
g2 = -5 + 3j
g3 = 2 + 7j
print(p_weierstrass_from_g2_g3(g2, g3)(z))
#(-1.36017042446944 - 2.49747290104859j)

def p_weierstrass_from_w1_w2(w1, w2):
    if w2.imag*w1.real < w1.imag*w2.real:
        raise ValueError("No solution. Do you want to exchange `w1 and `w2`?")
    ratio = w2 / w1
    if ratio.imag <= 0:
        raise ValueError("The ratio `w2/w1` must have a positive imaginary part.")
    q = exp(1j * pi * ratio)
    j2 = jtheta(2, 0, q)
    j3 = jtheta(3, 0, q)
    g2 = 4/3 * (pi/2/w1)**4 * (j2**8 - (j2*j3)**4 + j3**8) 
    g3 = 8/27 * (pi/2/w1)**6 * (j2**12 - (
        (3/2 * j2**8 * j3**4) + (3/2 * j2**4 * j3**8) 
    ) + j3**12)
    print("*******")
    print((g2, g3)) 
    print("*******")
    return p_weierstrass_from_g2_g3(g2, g3)

w1 = 1.0 + 4.0j
w2 = 1.0 + 15.0j
print(p_weierstrass_from_w1_w2(w1, w2)(z))
#(-0.0111828097206321 - 0.50738107410674j)

def p_weierstrass_from_tau(tau):
    return p_weierstrass_from_w1_w2(1.0, tau)

print(p_weierstrass_from_tau(1 + 4j)(z))
#(-0.430565782630798 - 3.62705469588323e-16j)

[stla]

And here is the Weierstrass Zeta function:

def zeta_weierstrass(g2, g3):
    r1, r2, r3 = polyroots([4, 0, -g2, -g3])
    a = sqrt(r1 - r3)
    b = sqrt(r1 - r2)
    c = sqrt(r2 - r3)
    if abs(a + b) < abs(a - b):
        b *= 1j
    if abs(a + c) < abs(a - c):
        c *= 1j
    if abs(c + 1j*b) < abs(c - 1j*b):
        a *= 1j 
        tmp = b
        b = 1j * c
        c = 1j * tmp
    w1 = 1 / agm(2*a, 2*b) * pi
    w3 = 1j / agm(2*a, 2*c) * pi
    print("#######################################")
    print(p_weierstrass_from_g2_g3(g2, g3)(w3)) # one of the ri
    q = exp(1j * pi * w3/w1)    
    p = 1 / w1 / 2
    eta1 = p / 6 / w1 * jtheta(1, 0, q, 3) / jtheta(1, 0, q, 1)
    return lambda z: (
        - eta1 * z * pi**2
        + pi*p * jtheta(1, pi*p*z, q, 1) / jtheta(1, pi*p*z, q)
    )

print("@@@@@@@@@@@@@@@@@@@@@@@@@@@@")
z = 1+1j
g2 = 5 + 3j
g3 = 5 + 3j
print(zeta_weierstrass(g2, g3)(z)) # same as Wolfram and R

[stla]

Here is a Costa surface. It is made with these two Weierstrass functions.

import pyvista as pv
import numpy as np

e1 = p_weierstrass_from_w1_w2(1/2, 1j/2)(1/2).real
c = 4*e1**2
zeta = np.frompyfunc(zeta_weierstrass(c, 0), 1, 1)
weier = np.frompyfunc(p_weierstrass_from_g2_g3(c, 0), 1, 1)
def fx(u, v):
    z = u + 1j*v
    zout = pi*(u + pi/4/e1) - zeta(z) + pi/2/e1*(zeta(z-1/2) - zeta(z-1j/2))
    return zout.real
fxufunc = np.frompyfunc(fx, 2, 1)
logufunc = np.frompyfunc(log, 1, 1)
def f(u, v):
    zx = fxufunc(u, v)
    zy = fxufunc(v, u)
    p = weier(u + 1j*v)
    return (
        zx/2,
        zy/2,
        sqrt(pi/2) * logufunc(abs((p-e1)/(p+e1)))/2
    )

tofloat = np.vectorize(float, otypes=["float"])
U, V = np.meshgrid(np.linspace(0.05, 0.95, 100), np.linspace(0.05, 0.95, 100))
x, y, z = f(U, V)
grid = pv.StructuredGrid(tofloat(x), tofloat(y), tofloat(z))
grid.plot(
    smooth_shading=True, color="darkred", specular=15    
)

[stla]

And here is the Weierstrass "P prime" function (derivative of the P function).

def pprime_weierstrass_from_g2_g3(g2, g3):
    r1, r2, r3 = polyroots([4, 0, -g2, -g3])
    a = sqrt(r1 - r3)
    b = sqrt(r1 - r2)
    c = sqrt(r2 - r3)
    w1 = None
    if abs(a + b) < abs(a - b):
        b *= -1
    if abs(a + c) < abs(a - c):
        c *= -1
    if abs(c + 1j*b) < abs(c - 1j*b):
        a = sqrt(r3 - r1)
        b = sqrt(r3 - r2)
        c = sqrt(r2 - r1)
        w1 = 1 / agm(1j*b, c) 
    else:
        w1 = 1 / agm(a, b)
    w3 = 1j / agm(a, c) 
    q = exp(1j * pi * w3/w1)
    f = jtheta(1, 0, q, 1)**3 / (jtheta(2, 0, q) * jtheta(3, 0, q) * jtheta(4, 0, q))
    s = 1 / w1 
    return lambda z: (
        -2*(1/w1)**3 * jtheta(2,s*z,q) * jtheta(3,s*z,q) * jtheta(4,s*z,q)
              * f / jtheta(1, s*z, q)**3
    )

z = 1+1j
g2 = 5 + 3j
g3 = 2 + 7j
print(pprime_weierstrass_from_g2_g3(g2, g3)(z))
# (-1.78040154378359 + 0.213468471404325j)

[stla]

Cleaner code here.

[stla]

As compared to Wikipedia, a factor pi is missing in your formula. I'm going to try.

[stla]

I understand now!! This equality is correct. Using my implementation, one gets it with $$\omega_1=1/2 \quad \text{and} \quad \omega_2=\tau/2.$$ I missed the one-half factor! Now I will try to understand how this generalizes to arbitrary omega's.

[stla]

Yes!! The relation is
$$\wp(z; ; \omega_1/2, ; \omega_2/2) = \wp(z/\omega_1; \tau) / \omega_1^2.$$
Your formula along with this relation highly simplify the implementation !!

[stla]

Almost finished.

You're quite silent :-) Are you still interested? Maybe you're busy.

I don't understand why the results I get when using the elliptical invariants (g2 and g3) are slightly different from the Wolfram results, e.g. 500.009741 vs 500.009715. With the R package elliptic, results are closer to the Wolfram results.

[stla]

That's fine, I found the error. Now I exactly get Wolfram!

[HeskethGD]

It would be great to get the Weierstrass elliptic functions into a package like mpmath and or sympy. Elliptic functions are applicable to may physical scenarios and yet seem rather neglected in Python packages. Trying to define them in terms of theta functions in every project would take ages. Great work @stla this helped me a lot!

[stla]

Hello @HeskethGD

I'm glad my code is helping someone. I did a package now (not published).

[HeskethGD]

Nice! If you're accepting feature requests the inverse of Weierstrass P would be awesome too! :) https://reference.wolfram.com/language/ref/InverseWeierstrassP.html

[stla]

I believe I implemented it in a R package. I'll port it to Python one day.