# Does Mathematica gives us a wrong result for the integral of a function including elliptic functions?

Calculus tells us that the differentiation of the integral of a function should be itself, but at least in one case Mathematica answers NO. I feel very confused. The new figure seems to bifurcate at one real root downward.

The story is that I try to integrate a function including a Jacobi elliptic function, then the strange thing happens. When I plot the figure of the integral, I found the slope sign does not agree with the sign of the integrand, that is, the original function. I cannot figure out what is wrong with Mathematica or if there is something tricky about elliptic functions and elliptic integrals? The original function and its integral by Mathematica software are shown as follows:

$$\text{mini}(m)= x_1 + \frac{x_2-x_1}{1-\frac{x_2-x_4}{x_1-x_4}\text{sn}^2\left[\frac{m-m_b}{\xi},k\right]}.$$ where $$\xi=\frac{4}{i\sqrt{\left(x_2-x_3\right)\left(x_1-x_4\right)}}$$ and the Jacobi modulus $$m=k^2=\frac{\left(x_1-x_3\right)\left(x_2-x_4\right)}{\left(x_2-x_3\right)\left(x_1-x_4\right)}$$, where $$x_1=-2.73205$$, $$x_2=0.732051$$, $$x_3=1 - i$$ and $$x_4=1+i$$. The constant $$m_b$$ can be calculated as $$m_b=\xi F\left[i \text{arcsinh}\left(\sqrt{-\frac{(0.5 - x_2)(x_1 - x_4)}{(0.5 - x_1)(x_2 - x_4)}}\right),k\right]=0.988254$$. Now we can plot the figure for the function $$\text{mini}(m)$$ with $$-10 as follows:

$$\int \text{mini}(m) \ \text{d}m = x_1 m - \frac{\xi\left(x_1-x_2\right)\Pi\left[n,\frac{m-m_b}{\xi},k\right]\text{dn}\left[\frac{m-m_b}{\xi},k\right]}{\sqrt{1-k^2\text{sn}^2\left[\frac{m-m_b}{\xi},k\right]}}.$$ where the elliptic characteristic $$n=\frac{x_2-x_4}{x_1-x_4}$$. Now we can plot the figure for the integral of the function $$\text{mini}(m)$$ with $$-10 as follows:

After calculating the numerical differentiation and plotting the figure by Mathematica, the figure shown below is not consistent with the original funcion $$\text{mini}(m)$$. But this third figure is consistent with the figure of the second one. It seems to show that Mathematica failed to calculate the correct integral of the original function $$\text{mini}(m)$$.

For the reference, my Mathematica code is attached as follows:

mini[m_] := (x1 x2 - x2 x4 - x1 x2 JacobiSN[1/4 I (m - mb) Sqrt[(x2 - x3) (x1 - x4)], ((x1 - x3) (x2 - x4))/((x2 - x3) (x1 - x4))]^2 + x1 x4 JacobiSN[
1/4 I (m - mb) Sqrt[(x2 - x3) (x1 - x4)], ((x1 - x3) (x2 - x4))/((x2 - x3) (x1 - x4))]^2)/((-x2 + x4) (-(x1/(x2 - x4)) + x4/(x2 - x4) + JacobiSN[1/4 I (m - mb) Sqrt[(x2 - x3) (x1 - x4)], ((x1 - x3) (x2 - x4))/((x2 - x3) (x1 - x4))]^2))
intmini[m_] := Integrate[mini[m], m]
intmini[m]
x1 = Root[#1^4 - 4*#1^2/1 + 8*#1/1 - 8*0.5 &, 1]
x2 = Root[#1^4 - 4*#1^2/1 + 8*#1/1 - 8*0.5 &, 2]
x3 = Root[#1^4 - 4*#1^2/1 + 8*#1/1 - 8*0.5 &, 3]
x4 = Root[#1^4 - 4*#1^2/1 + 8*#1/1 - 8*0.5 &, 4]
mb = -((4 I EllipticF[I*ArcSinh[Sqrt[-(((0.5 - x2)*(x1 - x4))/((0.5 - x1)*(x2 - x4)))]], ((x1 - x3)*(x2 - x4))/((x2 - x3)*(x1 - x4))])/Sqrt[(x2 - x3)*(x1 - x4)])
Plot[Re[mini[m]], {m, -10, 10}]
Plot[(m - mb) x1 + (4 I (x1 - x2) EllipticPi[(x2 - x4)/(x1 - x4), JacobiAmplitude[1/4 I (m - mb) Sqrt[(x2 - x3) (x1 - x4)], ((x1 - x3) (x2 - x4))/((x2 - x3) (x1 - x4))], ((x1 - x3) (x2 - x4))/((x2 - x3) (x1 - x4))] JacobiDN[1/4 I (m - mb) Sqrt[(x2 - x3) (x1 - x4)], ((x1 - x3) (x2 - x4))/((x2 - x3) (x1 - x4))])/(Sqrt[(x2 - x3) (x1 - x4)] Sqrt[1 + ((x1 - x3) (x2 - x4) JacobiSN[1/4 I (m - mb) Sqrt[(x2 - x3) (x1 - x4)], ((x1 - x3) (x2 - x4))/((x2 - x3) (x1 - x4))]^2)/((-x2 + x3) (x1 - x4))]), {m, -10, 10}]
Needs["NumericalCalculus"]
Plot[Re[ND[(m - mb) x1 + (4 I (x1 - x2) EllipticPi[(x2 - x4)/(x1 - x4), JacobiAmplitude[1/4 I (m - mb) Sqrt[(x2 - x3) (x1 - x4)], ((x1 - x3) (x2 - x4))/((x2 - x3) (x1 - x4))], ((x1 - x3) (x2 - x4))/((x2 - x3) (x1 - x4))] JacobiDN[1/4 I (m - mb) Sqrt[(x2 - x3) (x1 - x4)], ((x1 - x3) (x2 - x4))/((x2 - x3) (x1 - x4))])/(Sqrt[(x2 - x3) (x1 - x4)] Sqrt[1 + ((x1 - x3) (x2 - x4) JacobiSN[1/4 I (m - mb) Sqrt[(x2 - x3) (x1 - x4)], ((x1 - x3) (x2 - x4))/((x2 - x3) (x1 - x4))]^2)/((-x2 + x3) (x1 - x4))]), m, p]], {p, -10, 10}]


Again, I would like to comment that the first plot is for the figure of the original function, that is, the integrand $$\text{mini}(m)$$. The second the plot is for the figure of its integral. The third one is for the figure of the differentiation of the integral. The third figure is expected to be the same as the first one, but it is different. The signs in the third figure agrees with the signs in the second figure of the integral. It seems to tell us that Mathematica gives a wrong integral for the function including elliptic functions, like the original function I used. I am not sure how to calculate a correct integral of my original function. Welcome any useful suggestions and helps in Math and Mathematica!

PS: I would like to add one additional comment. In my case here, the elliptic modulus m and Jacobi amplitude are both complex numbers, not real numbers. I am wondering if there is any general theory of elliptic functions and elliptic integrals beyond the real elliptic modulus and real Jacobi amplitude?

• No, I think ND[ ] is correct, because Fig 3 in the post is consistent with Fig 2 in the post. But the problem is that the plot of the integral (Fig 2) is not consistent with the plot of the integrand (Fig 1). There should be more oscillations along the curve in Fig 2. It seems to show that Mathematica obtains a wrong integral of the function mini(x) – Hao Wu Jan 16 '19 at 3:37

Working with exact arithmetic (until the end) indicates this example at least works out as expected.

mini = (x1 x2 - x2 x4 -
x1 x2 JacobiSN[
1/4 I (m -
mb) Sqrt[(x2 - x3) (x1 - x4)], ((x1 - x3) (x2 - x4))/((x2 -
x3) (x1 - x4))]^2 +
x1 x4 JacobiSN[
1/4 I (m -
mb) Sqrt[(x2 - x3) (x1 - x4)], ((x1 - x3) (x2 - x4))/((x2 -
x3) (x1 - x4))]^2)/((-x2 + x4) (-(x1/(x2 - x4)) +
x4/(x2 - x4) +
JacobiSN[
1/4 I (m -
mb) Sqrt[(x2 - x3) (x1 - x4)], ((x1 - x3) (x2 -
x4))/((x2 - x3) (x1 - x4))]^2));

intmini = Integrate[mini, m];
diff = Simplify[D[intmini, m] - mini];


rep1 = {x1 -> Root[#1^4 - 4*#1^2/1 + 8*#1/1 - 4 &, 1],
x2 -> Root[#1^4 - 4*#1^2/1 + 8*#1/1 - 4 &, 2],
x3 -> Root[#1^4 - 4*#1^2/1 + 8*#1/1 - 4 &, 3],
x4 -> Root[#1^4 - 4*#1^2/1 + 8*#1/1 - 4 &, 4]};
rep2 = {mb ->
Simplify[-((4 I EllipticF[
I*ArcSinh[
Sqrt[-(((1/2 - x2)*(x1 - x4))/((1/2 - x1)*(x2 -
x4)))]], ((x1 - x3)*(x2 - x4))/((x2 - x3)*(x1 -
x4))])/Sqrt[(x2 - x3)*(x1 - x4)]) /. rep1]};


Use them to reduce that difference to a univariate function.

mdiff = Simplify[diff /. rep1 /. rep2]

(* Out[4797]= -1 - Sqrt[3] + (Sqrt[
3] ((-2 - I) + Sqrt[3]) ((2 - I) + Sqrt[3]) EllipticPi[((2 + I) -
Sqrt[3])/((2 + I) + Sqrt[3]),
JacobiAmplitude[(I m)/(I + Sqrt[3]) -
EllipticF[
I ArcSinh[Sqrt[(-3 I + (1 + 2 I) Sqrt[3])/(
3 I + (1 + 2 I) Sqrt[3])]], (I - Sqrt[3])/(I + Sqrt[3])], (
I - Sqrt[3])/(I + Sqrt[3])], (I - Sqrt[3])/(
I + Sqrt[3])] JacobiCN[(I m)/(I + Sqrt[3]) -
EllipticF[
I ArcSinh[Sqrt[(-3 I + (1 + 2 I) Sqrt[3])/(
3 I + (1 + 2 I) Sqrt[3])]], (I - Sqrt[3])/(I + Sqrt[3])], (
I - Sqrt[3])/(
I + Sqrt[3])] JacobiDN[(I m)/(I + Sqrt[3]) -
EllipticF[
I ArcSinh[Sqrt[(-3 I + (1 + 2 I) Sqrt[3])/(
3 I + (1 + 2 I) Sqrt[3])]], (I - Sqrt[3])/(I + Sqrt[3])], (
I - Sqrt[3])/(
I + Sqrt[3])]^2 JacobiSN[(I m)/(I + Sqrt[3]) -
EllipticF[
I ArcSinh[Sqrt[(-3 I + (1 + 2 I) Sqrt[3])/(
3 I + (1 + 2 I) Sqrt[3])]], (I - Sqrt[3])/(I + Sqrt[3])], (
I - Sqrt[3])/(I + Sqrt[3])])/((-1 + I Sqrt[3]) ((1/(
I + Sqrt[
3]))(I + Sqrt[
3] + (-I + Sqrt[3]) JacobiSN[(I m)/(I + Sqrt[3]) -

EllipticF[
I ArcSinh[Sqrt[(-3 I + (1 + 2 I) Sqrt[3])/(
3 I + (1 + 2 I) Sqrt[3])]], (I - Sqrt[3])/(
I + Sqrt[3])], (I - Sqrt[3])/(I + Sqrt[3])]^2))^(
3/2)) + (Sqrt[
3] ((1 + 2 I) + I Sqrt[3]) ((-2 - I) + Sqrt[
3]) EllipticPi[((2 + I) - Sqrt[3])/((2 + I) + Sqrt[3]),
JacobiAmplitude[(I m)/(I + Sqrt[3]) -
EllipticF[
I ArcSinh[Sqrt[(-3 I + (1 + 2 I) Sqrt[3])/(
3 I + (1 + 2 I) Sqrt[3])]], (I - Sqrt[3])/(I + Sqrt[3])], (
I - Sqrt[3])/(I + Sqrt[3])], (I - Sqrt[3])/(
I + Sqrt[3])] JacobiCN[(I m)/(I + Sqrt[3]) -
EllipticF[
I ArcSinh[Sqrt[(-3 I + (1 + 2 I) Sqrt[3])/(
3 I + (1 + 2 I) Sqrt[3])]], (I - Sqrt[3])/(I + Sqrt[3])], (
I - Sqrt[3])/(
I + Sqrt[3])] JacobiSN[(I m)/(I + Sqrt[3]) -
EllipticF[
I ArcSinh[Sqrt[(-3 I + (1 + 2 I) Sqrt[3])/(
3 I + (1 + 2 I) Sqrt[3])]], (I - Sqrt[3])/(I + Sqrt[3])], (
I - Sqrt[3])/(
I + Sqrt[3])])/((I + Sqrt[3]) \[Sqrt]((1/(
I + Sqrt[
3]))(I + Sqrt[
3] + (-I + Sqrt[3]) JacobiSN[(I m)/(I + Sqrt[3]) -

EllipticF[
I ArcSinh[Sqrt[(-3 I + (1 + 2 I) Sqrt[3])/(
3 I + (1 + 2 I) Sqrt[3])]], (I - Sqrt[3])/(
I + Sqrt[3])], (I - Sqrt[3])/(I + Sqrt[3])]^2))) - ((1 +
I) (-I + Sqrt[
3] + (I + Sqrt[3]) JacobiSN[(I m)/(I + Sqrt[3]) -
EllipticF[
I ArcSinh[Sqrt[(-3 I + (1 + 2 I) Sqrt[3])/(
3 I + (1 + 2 I) Sqrt[3])]], (I - Sqrt[3])/(
I + Sqrt[3])], (I - Sqrt[3])/(I + Sqrt[3])]^2))/((2 + I) +
Sqrt[3] + ((-2 - I) + Sqrt[3]) JacobiSN[(I m)/(I + Sqrt[3]) -
EllipticF[
I ArcSinh[Sqrt[(-3 I + (1 + 2 I) Sqrt[3])/(
3 I + (1 + 2 I) Sqrt[3])]], (I - Sqrt[3])/(I + Sqrt[3])], (
I - Sqrt[3])/(I + Sqrt[3])]^2) + ((1/4 + I/4) Sqrt[
3] ((-2 + I) + Sqrt[3]) ((2 + I) + Sqrt[3])^2 JacobiDN[(I m)/(
I + Sqrt[3]) -
EllipticF[
I ArcSinh[Sqrt[(-3 I + (1 + 2 I) Sqrt[3])/(
3 I + (1 + 2 I) Sqrt[3])]], (I - Sqrt[3])/(I + Sqrt[3])], (
I - Sqrt[3])/(I + Sqrt[3])]^2)/(-1 - Sqrt[
3] - (1 - 3 I) JacobiSN[(I m)/(I + Sqrt[3]) -
EllipticF[
I ArcSinh[Sqrt[(-3 I + (1 + 2 I) Sqrt[3])/(
3 I + (1 + 2 I) Sqrt[3])]], (I - Sqrt[3])/(I + Sqrt[3])], (
I - Sqrt[3])/(
I + Sqrt[3])]^2 + (-1 + Sqrt[3]) JacobiSN[(I m)/(I + Sqrt[3]) -
EllipticF[
I ArcSinh[Sqrt[(-3 I + (1 + 2 I) Sqrt[3])/(
3 I + (1 + 2 I) Sqrt[3])]], (I - Sqrt[3])/(I + Sqrt[3])], (
I - Sqrt[3])/(I + Sqrt[3])]^4) *)


Evaluate that difference on the path used for those plots:

Max[Abs[N[Table[mdiff, {m, -10, 10, 1/10}], 20]]]

(* During evaluation of In[4801]:= N::meprec: Internal precision limit \$MaxExtraPrecision = 50. reached while evaluating {-1-Sqrt[3]-<<1>>/((I+Sqrt[3]) Sqrt[<<1>>/(3 <<1>>)])-(Sqrt[3] <<5>> <<1>>)/((-1+I <<1>>) <<1>>)-((1+I) (-I+Sqrt[3]+(I+Sqrt[3]) JacobiSN[<<1>>,<<1>>]^2))/((2+I)+Sqrt[3]+((-2-I)+Power[<<2>>]) JacobiSN[<<2>>]^2)+((1/4+I/4) Sqrt[3] ((-2+I)+Sqrt[3]) ((2+I)+Sqrt[3])^2 JacobiDN[10 I Power[<<2>>]+<<1>>,(I<<1>>Times<<1>><<1>>])/Plus[<<2>>]]^2)/(-1-Sqrt[3]-(1-3 I) JacobiSN[<<2>>]^2+(-1+Power[<<2>>]) JacobiSN[<<2>>]^4),<<1>>,<<47>>,<<1>>,<<151>>}.

Out[4801]= 0.*10^-56 *)


So the evaluations are, as best one can tell, all zero.

• The problem is that the plot of the integral (Fig 2 in the post) is not consistent with the plot of the integrand (Fig 1 in the post). There should be more oscillations along the curve in Fig 2. It seems to show that Mathematica obtains a wrong integral of the function mini(x) – Hao Wu Jan 16 '19 at 3:34

This is just a guess, but the problem seems to be that$$\sqrt{x^2}\ne x$$ in general. The integral according to Mathematica had to choose a square root and it is wrong half of the time. More precisely, the integral that Mathematica returns has in the denominator a factor something like $$\,\sqrt{1 - m\, \textrm{sn}(x,m)^2}\,$$ and$$\,\textrm{dn}(x,m)\,$$ in the numerator. The two functions are equal but sometimes the wrong square root is chosen. After this is fixed, the results are correct. But I could be wrong.

Here is my general code:

ClearAll[x1, x2, x3, x4, m, ma, mb, mc, md, mx, sna2,
mini, imini, iminigood, dimini, dimingood, eps];
ma = ((x1 - x3) (x2 - x4))/((x2 - x3) (x1 - x4));
mc = (x2 - x4)/(x1 - x4);
md = Sqrt[(x2 - x3) (x1 - x4)]; mx = 1/4 I md;

sna2 = JacobiSN[(m - mb) mx, ma]^2;
mini = (x2 - mc*x1*sna2)/(1 - mc*sna2);
imini = Integrate[mini, m];
iminigood = x1*(m - mb) - (x1 - x2)/mx *
EllipticPi[mc, JacobiAmplitude[(m - mb) mx, ma], ma];
dimini = ((imini /. m -> m + eps) - imini)/eps;
diminigood = ((iminigood /. m -> m + eps) - iminigood)/eps;


and for specific numerical values:

x1 = Root[#1^4 - 4*#1^2/1 + 8*#1/1 - 8*0.5 &, 1];
x2 = Root[#1^4 - 4*#1^2/1 + 8*#1/1 - 8*0.5 &, 2];
x3 = Root[#1^4 - 4*#1^2/1 + 8*#1/1 - 8*0.5 &, 3];
x4 = Root[#1^4 - 4*#1^2/1 + 8*#1/1 - 8*0.5 &, 4];
mb = Re[-((4 I EllipticF[ I*ArcSinh[
Sqrt[-(((0.5 - x2)*(x1 - x4))/((0.5 - x1)*
(x2 - x4)))]], ((x1 - x3)*(x2 - x4))/((x2 - x3)*
(x1 - x4))])/Sqrt[(x2 - x3)*(x1 - x4)])];
eps = .000001;
Print[{ma, md} - {Exp[2 Pi I/3], Sqrt[3] - I} // N // Chop];
{ma, md} = {Exp[2 Pi I/3], Sqrt[3] - I} // N;


Of course, the iminigood has jump discontinuities, but the fix is to use

iminigood + Quotient[m + c1, c2] c3


with approximate values of c1 = 2.3832, c2 = 6.7429, c3 = 15.511

Here is a plot with this fix added:

Plot[iminigood + Quotient[m + 2.3832, 6.7429] 15.511, {m, -10, 10}]


More work could be done to pin down why the jump discontinuities appear. Clearly, the integral of a smooth function is smooth. Something about the EllipticPi[] is not quite right, but I don't yet know what it is.

• Thanks! This is helpful information. But I have a further question: why the integral of a smooth function can become not smooth? For example, after integrating the mini(x), the incomplete elliptic integral of the third kind shows up. Its plot is not smooth. Does it mean that the integral can change the smoothness of the integrand? It can make the curve not smooth. When I cancel dn and sqrt(1-m sn^2) with the assumption of the same sign or oppose sign, only the incomplete elliptic integral of the third kind is left. And then plot it, I find that the curve is jumping periodically, like tan(x). – Hao Wu Jan 16 '19 at 21:03
• Thanks! Yes, I also wonder EllipticPi[ ] in Mathematica seems not right and gives the discontinuities in the plot. – Hao Wu Jan 21 '19 at 4:15
• what is "imini2" in the second part of your code? – Hao Wu Jan 21 '19 at 15:33
• And how did you fix it? I am wondering if this way will change the original slope? – Hao Wu Jan 21 '19 at 15:36
• "iminigood" was "imini2" before I renamed it. Fixed now. – Somos Jan 21 '19 at 15:47