# How to evaluate theta function's derivative numerically?

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...

• Which to accept will depend on the vote, for me both of them are excellent. – an offer can't refuse Nov 2 '19 at 7:16

### 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  *)

• I fail to see how your formula gives the derivative. In the Cauchy formula for the derivative the $z-z_0$ should be squared, and there should be $2\pi i$ instead of $2\pi$. What do I miss? – mickep Nov 2 '19 at 12:40
• @mickep I forgot the square and the $i$ -- thanks! (And the $dz$.) – Michael E2 Nov 2 '19 at 13:54
• Thanks, and nice! I was confused there for a while! – mickep Nov 2 '19 at 15:35
• As it turns out, there is an undocumented function for expanding out the genus-1 case of SiegelTheta[] in terms of EllipticTheta[]: D[SiegelTheta[{{a}, {b}}, {{t}}, z] // SystemSiegelThetaDumpSimplifySiegelTheta, z]. – J. M.'s discontentment Mar 2 at 5:43

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

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
`