How to evaluate theta function's derivative numerically?

Question

I ran into this derivative that Mathematica won't evaluate:

ClearAll[f, g];
a = 0.1;
b = 0.2;
t = Exp[I 2 Pi/3];
f[z_] := SiegelTheta[{{a}, {b}}, {{t}}, z]
g[z_] = D[f[z], z]
g[0.1] // N

I am not sure if this derivative can be evaluated to a numerical value actually...

Answer 1

Numerical derivative

Based on $$f'(z_0)={1 \over 2\pi i}\,\int_\gamma {f(z) \, dz\over (z-z_0)^2}\,,$$ where $\gamma$ is a closed contour containing $z_0$ in its interior.

fPrime[z0_] :=
  1/8 Sum[f[z0 + dz]/dz,
          {dz, Exp[2 Pi I Most@Subdivide[0., 1., 8]]/1000}];

fPrime[0.1]
(*  -0.256724 + 1.47096 I  *)

---

Update:

Discretizing the integral with n = 2 points instead of n = 8 yields the central difference formula, and for a radius of Abs[dz] == 1*^-9, it will have a truncation error less than machine-precision for analytic functions whose higher-order derivatives do not grow too rapidly. To prevent round-off error overwhelming the truncation error, we compute f[z] at high precision. This is faster than the 8-point machine-precision code above on the OP's function (I suspect because SiegelTheta is somewhat expensive to compute). The 8-point formula with a radius of 1/1000 in fPrime has a relative error of $10^{-10}$ or less in a neighborhood of $z = 0.1 + 0i$. The function ND[] has a relative error of $10^{-5}$ or less. Over the square with ReIm[z] between ±1, the relative errors of fPrime and ND can be a couple of orders of magnitude larger, but fPrime2 below maintains machine-precision-accurate results.

ClearAll[f, fPrime2];
a = 1/10;
b = 2/10;
t = Exp[I 2 Pi/3];
f[z_] := SiegelTheta[{{a}, {b}}, {{t}}, z]
fPrime2[z_?NumericQ] := N@With[{z0 = SetPrecision[z, 32], r = 1*^-9},
    (f[r + z0] - f[-r + z0])/(2 r)
    ];

Symbolic derivative

For the OP's special case of SiegelTheta[], a symbolic derivative can be computed from the Sum[] of its theta series expansion, which returns a sum in terms of EllipticTheta[], whose derivative is implemented as EllipticThetaPrime[[]:

SiegelThetaPrime[{{a_}, {b_}}, {{t_}}, z_] = Simplify@D[
   Sum[Exp[
     I Pi ((n + {a}).{{t}}.(n + {a}) + 
        2 (n + {a}).(z + {b}))], {n, -Infinity, Infinity}],
   z]
(* 
(E^(-((I π (b + z)^2)/
  t)) π (-2 I (b + z) EllipticTheta[3, (π (b + a t + z))/t, 
     E^(-((I π)/t))] + 
   EllipticThetaPrime[3, (π (b + a t + z))/t, 
    E^(-((I π)/t))]))/(Sqrt[-I t] t)
*)

SiegelThetaPrime[{{1/10}, {1/5}}, {{Exp[I 2 π/3]}}, 0.1]
(*  -0.256724 + 1.47096 I  *)

Answer 2

You can compute a numerical derivative as follows:

ClearAll[f, g];
Needs["NumericalCalculus`"]
a = 0.1;
b = 0.2;
t = Exp[I 2 Pi/3];
f[z_] := SiegelTheta[{{a}, {b}}, {{t}}, z]
g[z0_] := ND[f[z], z, z0]
g[0.1] 

(*-0.256725 + 1.47096 I*)

I haven't checked the result is correct.

Answer 3

A quick-and dirty method is to use complex-step differentiation:

With[{a = 1/10, b = 1/5, t = Exp[I 2 π/3], z = 1/10, h = 10^-9, prec = 20}, 
     N[(SiegelTheta[{{a}, {b}}, {{t}}, z + I h] -
        SiegelTheta[{{a}, {b}}, {{t}}, z - I h])/(2 I h), prec]]
   -0.2567239264794337275 + 1.4709617732598025465 I

where even a modest-sized step size can yield a slightly more accurate result, compared to using a purely real step size.

Alternatively, one can use Cauchy's differentiation formula. Michael's answer shows one possible implementation, and here is another one:

With[{a = 1/10, b = 1/5, t = Exp[I 2 π/3], z = 1/10, r = 10^-6}, 
     NIntegrate[SiegelTheta[{{a}, {b}}, {{t}}, z + r Exp[I u]]/(2 π r Exp[I u]),
                {u, -π, π}, Method -> "Trapezoidal", WorkingPrecision -> 20]]
   -0.2567239264794337266 + 1.4709617732598025411 I

Finally, one might also consider trying the "Lanczos derivative":

With[{a = 1/10, b = 1/5, t = Exp[I 2 Pi/3], z = 1/10, h = 10^-9},
     (3/(2 h^3)) NIntegrate[u SiegelTheta[{{a}, {b}}, {{t}}, z + u], {u, -h, h}, 
                            Method -> "GlobalAdaptive", WorkingPrecision -> 20]]
   -0.25672392647943372659 + 1.4709617732598025411 I