# Symbolic integration of elliptic functions

Is there a clever way to integrate products of elliptic functions like $\wp(z;g_2,g_3)$ or $\zeta(z;g_2,g_3)$ in Mathematica?

Mathematica seems to be able to integrate functions like

Integrate[WeierstrassP[z, {g2, g3}]^3, z]


However, it cannot integrate

Integrate[WeierstrassP[z, {g2, g3}] WeierstrassP[z + a, {g2, g3}], z]


I have tried using the addition formulae for $\wp(z)$, but Mathematica is still unable to give the answer.

Unfortunately, as I have said every so often, Mathematica's symbolic handling of elliptic functions is far from optimal. To evaluate the OP's second integral, we will need some heavy machinery.

N. I. Akhiezer, in his book, reveals a formula that decomposes an arbitrary elliptic function into a so-called "partial fraction" form that exposes the poles of the elliptic function. Mathematica's usefulness here is in its ability to expand the Weierstrass functions in series, which we can use to obtain the desired form.

f[z_] := WeierstrassP[z, {g2, g3}] WeierstrassP[z + a, {g2, g3}]


At once, it can be shown that f[z] has poles at $z=0$ and $z=-a$. Consideration of the function's double periodicity shows that all the other poles of f[z] are congruent to these two, so it suffices to consider only these two in the "partial fraction" form.

To derive the formula, we start by taking the "principal part" of the series expansion of f[z] at the two poles:

Series[f[z], {z, 0, -1}] // Normal
WeierstrassP[a, {g2, g3}]/z^2 + WeierstrassPPrime[a, {g2, g3}]/z

Series[f[z], {z, -a, -1}] // Normal
WeierstrassP[a, {g2, g3}]/(a + z)^2 - WeierstrassPPrime[a, {g2, g3}]/(a + z)


We can now construct the tentative "partial fraction" representation like so (compare the expression here with the one in Akhiezer's book):

ft[z_] := WeierstrassPPrime[a, {g2, g3}] WeierstrassZeta[z, {g2, g3}] +
WeierstrassP[a, {g2, g3}] WeierstrassP[z, {g2, g3}] +
(-WeierstrassPPrime[a, {g2, g3}]) WeierstrassZeta[z + a, {g2, g3}] +
WeierstrassP[a, {g2, g3}] WeierstrassP[z + a, {g2, g3}]


As noted by Akhiezer, f[z] and ft[z] differ by a constant. To determine this constant, we use Series[] again:

c = Normal[Series[f[z] - ft[z], {z, 0, 0}]]
-g2/4 + 2 WeierstrassP[a, {g2, g3}]^2 +
WeierstrassPPrime[a, {g2, g3}] WeierstrassZeta[a, {g2, g3}]


The final "partial fraction" decomposition of f[z] is thus c + ft[z].

Going back to the OP's problem, the reason for converting to the "partial fraction" form is similar to the reason for converting arbitrary rational functions to partial fractions: it eases the task of integration.

Collect[Integrate[c + ft[z], z], z]
Log[WeierstrassSigma[z, {g2, g3}]] WeierstrassPPrime[a, {g2, g3}] -
Log[WeierstrassSigma[a + z, {g2, g3}]] WeierstrassPPrime[a, {g2, g3}] +
z (-g2/4 + 2 WeierstrassP[a, {g2, g3}]^2 +
WeierstrassPPrime[a, {g2, g3}] WeierstrassZeta[a, {g2, g3}]) -
WeierstrassP[a, {g2, g3}] WeierstrassZeta[z, {g2, g3}] -
WeierstrassP[a, {g2, g3}] WeierstrassZeta[a + z, {g2, g3}]


One could probably simplify this result further by using the addition formulae for the Weierstrass functions, but I'll leave that to you.