I am trying to calculate/plot the derivative of the second Jacobi theta function $d\theta_2(0, e^{-\pi t} )/dt$.

Calculating or plotting the function itself works fine:

In[1]:= theta2[t_] := EllipticTheta[2, 0, Exp[-1*Pi*t]];

Out[2]= 1.07398

with the result being real - as expected from Mathematica's definition:

enter image description here

However, when I try to calculate derivatives of the above, I get a significant imaginary part:

In[3]:= dtheta2[t_] = D[EllipticTheta[2, 0, Exp[-1*Pi^2 *t]], t];

Out[4]= -0.794774 + 0.280078 I

(Using Set (=) rather than SetDelayed (:=) as discussed. Taking the derivative of theta2[t] instead does not seem to make a difference).

Any ideas what might be going on?

(Note, the branch cut is taken from 0 to -1, so should not be an issue.)

  • $\begingroup$ Welcome to Mathematica.SE! I suggest the following: 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Read the faq! 3) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! $\endgroup$ – user9660 Apr 21 '15 at 5:03

The complex results seem to stem from a peculiarity of Mathematica's implementation of EllipticTheta[] and/or Derivative[]. To demonstrate this, let's define the derivative with respect to the third argument:

ϑ2p[q_] := Derivative[0, 0, 1][EllipticTheta][2, 0, q]

Now an innocent question: what is the numerical value of ϑ2p[1/20]? Let's try:

   6.202438815596758 - 2.6108233971330863 I

It's complex?! Really? But what if we use arbitrary precision?

N[ϑ2p[1/20], 20]
   1.3160596143868681205 + 2.5399913229586444864 I

N[ϑ2p[1/20], 25]
   1.316059614386868120466950 + 2.539991322958644486424198 I

N[ϑ2p[1/20], 30]
   0.84478842638869011647966406334 - 2.86323146018455201167660171753 I

N[ϑ2p[1/20], 40]
   3.472142008910041601254984400930093009741 -
   2.927871724716208424440402760713828235294 I

By Jove! Where, indeed, has the consistency gone? Now, which of these results are correct?

To see if we can regain some of our sanity back (doubtful, but let's try anyway), we'll use the classical symmetric difference as a verification:

ϑdiff[q_, h_] := (EllipticTheta[2, 0, q + h] - EllipticTheta[2, 0, q - h])/(2 h)

Table[N[ϑdiff[1/20, 2^-k], 20], {k, 6}]
   {0.8898189458073228447 - 1.4025445411536147816 I,
    1.2616570729887737889 - 1.9672636094000931928 I,
    2.3559020961062510658 - 2.9769836079456671254 I,
    5.6000758270930682868 - 3.7835575245327611098 I,
    5.3519419508807769750, 4.9421017280379492149}

Hmm, the first few values are complex. Since Derivative[] is known to do high-order differencing in numerical evaluations, this may be why we're getting complex results. Nevertheless, the last few values are real (since 1/20 - 2^-5 and the terms beyond are positive), so let's continue the series further:

diffs = Table[N[ϑdiff[1/20, 2^-k], 20], {k, 5, 12}]
   {5.3519419508807769750, 4.9421017280379492149, 4.8608185742381784551,
    4.8414731561104121633, 4.8366939097686808577, 4.8355026115339014653,
    4.8352050057014693893, 4.8351306179004033988}

At a higher precision, I now get consistent results, so I am now somewhat confident that this sequence has the semblance of convergence. To get a better estimate, I shall now attempt Richardson extrapolation:

InterpolatingPolynomial[Transpose[{2^-Range[5, 12], diffs}], 0]

What if we go further, and use higher precision along the way (to compensate for the inevitable cancellation of digits inherent in the method)?

InterpolatingPolynomial[Table[{2^-k, N[ϑdiff[1/20, 2^-k], 30]}, {k, 5, 20}], 0]

InterpolatingPolynomial[Table[{2^-k, N[ϑdiff[1/20, 2^-k], 40]}, {k, 5, 28}], 0]

Ahh... much better. Unfortunately, this was all rather ad-hoc, so I don't know if you can bundle this strategy into a nice tidy routine whose guts you can forget about after writing it. But, it was a fun little experiment.


You are using insufficient precision:


(* -0.794774 + 0.280078 I *)


(* MachinePrecision *)

Use exact arguments:


and get the exact result. Use N et. al. with desired precision to retrieve numeric values.

N[theta2'[8/10], 10]

(* -0.8874928427 + 4.596*10^-7 I *)
  • $\begingroup$ Thx for the reply, @rasher. I've tried increasing precision before and did not get any substantive improvements wrt size of the imaginary term (e.g., with above I get: dtheta2[t_] = D[EllipticTheta[2, 0, Exp[-1*Pi *t]], t]; N[dtheta2[8/10], 10] N[theta2'[8/10], 10] output: -0.3317952035 - 0.4122258903 I -0.3317952035 - 0.4122258903 I) though the values in this case are different, which is also troubling...(using Mathematica 9). $\endgroup$ – albert818 Apr 21 '15 at 5:39
  • $\begingroup$ @albert818: Hmm. Something fishy perhaps - I restarted kernel, got a different result (not matching mine in post, or yours). Then did exact calculation, then applied N to that, got a different result, but when re-doing the N[...] surrounding all version, it now matches the latter, though both are different from the post and yours. It appears to be caching results but something gets whacked (try doing it with N[...,10], then re-do with say N[...,100], then re-do with N[...,10] again. Weird. Also on 9, btw. $\endgroup$ – ciao Apr 21 '15 at 6:00
  • $\begingroup$ Same story: did a complete program restart and got: -1.2947074204 - 0.4939796009 I with N[...,10], then ran N[...,100] and got -1.0078730012 + 0.0790778685 I (truncated paste). Back to N[...,10] - now same as N[...,100] (must be caching). Finally, did N[...,1000] and got yet another number and now same at N[...,10]: -0.8874927502 - 2.184*10^-7 I. This one has the Im part at close to chopping levels, at last(!) - but not sure I can believe it due to earlier instabilities. $\endgroup$ – albert818 Apr 21 '15 at 6:26
  • $\begingroup$ Tried plotting with N[...,500] while setting $MaxExtraPrecision = 500, but still get the same pattern of significant oscillatory Im component (though the frequency of oscillations seems to have gone up, which suggests numerics to me). The whole thing is particularly odd, since I believe theta series are generally rather fast-converging due to the n^2 term in the exponent... Any additional clues would be appreciated. $\endgroup$ – albert818 Apr 21 '15 at 7:35
  • $\begingroup$ @albert818: I'm at a loss to further theories - hopefully one of the wolfram gang that hangs out here might have some insight, and perhaps pinging wolfram support might get some info. $\endgroup$ – ciao Apr 21 '15 at 19:27

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.