5
$\begingroup$

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

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

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  *)
| improve this answer | |
$\endgroup$
  • $\begingroup$ 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? $\endgroup$ – mickep Nov 2 '19 at 12:40
  • $\begingroup$ @mickep I forgot the square and the $i$ -- thanks! (And the $dz$.) $\endgroup$ – Michael E2 Nov 2 '19 at 13:54
  • $\begingroup$ Thanks, and nice! I was confused there for a while! $\endgroup$ – mickep Nov 2 '19 at 15:35
  • $\begingroup$ Indeed, one can express the univariate case of SiegelTheta[] in terms of the more conventional Jacobi theta functions. $\endgroup$ – J. M.'s discontentment Nov 13 '19 at 13:23
  • 1
    $\begingroup$ 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] // System`SiegelThetaDump`SimplifySiegelTheta, z]. $\endgroup$ – J. M.'s discontentment Mar 2 at 5:43
6
$\begingroup$

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

| improve this answer | |
$\endgroup$
1
$\begingroup$

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
| improve this answer | |
$\endgroup$

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.