# Directional derivative of SiegelTheta

I'm working on a problem where I have to integrate both the Mathematica function SiegelTheta and some of its second order directional derivatives. Using the function works well but something goes wrong when I try to differentiate it. What I want to do is to take the second order derivative with respect to the vector $s$ in SiegelTheta[$\Omega$,s] and then dot it twice with the vectors (1,0),(1,0) or (0,1),(0,1). I want $s$ to be a two-component vector and $\Omega$ to be the matrix

$$\Omega(x,y,z) = \frac{i\sqrt{1/z - 1}}{\pi}\begin{pmatrix} (1 - y) & (-1+x+y) \\ (-1+x+y) & (1-x) \end{pmatrix}.$$

and after differentiating I want to let $s\to0$. Since I don't know how to take the derivative with respect to a vector I've tried doing it component-wise, using

S[Omega] = D[SiegelTheta[Omega,{s1,0}],{s1,2}]/.t1->0


and similarly for the derivative with respect to the second component of $s$. My problem is that I can't get the resulting function to take numerical values, even when evaluating it in a single point such as $(0.5,0.5,0.5)$. Upon integrating it I also get an error message saying the integrand takes non-numerical values in the whole region $x\in[0,1],y\in[0,1-x],z\in[0,1]$. To be more specific regarding the integration what I've tried using is

NIntegrate[S[Omega[x,y,z]],{x,0,1},{y,0,1-x},{z,0,1}].


My question is thus this: Do I have to somehow specify before taking the derivative that I want $\Omega$ to be a $2\times 2$ matrix and $s$ a vector with two components, and in that case how do I do that? Also, is there something happening during the differentiation that I'm not aware of that stops the function from taking numerical value? And is there something wrong with my approach to the integration?

I'm far from an expert on Mathematica and would appreciate any help I could get.

• Can you post the Mathematica expression for Ω – ssch Feb 14 '13 at 10:22
• I've used the expression $\Omega$ = {{(1/z - 1)^(1/2)*I/Pi*(1 - y), (1/z - 1)^(1/2)* I/Pi*(-1 + x + y)}, {(1/z - 1)^(1/2)* I/Pi*(-1 + x + y), (1/z - 1)^(1/2)*I/Pi*(1 - x)}} – Emil Bostrom Feb 14 '13 at 10:52
• What's the domain for Ω? When I try things like SiegelTheta[Ω[1, 3, 3], {1, 0, 1}] I just get the error SiegelTheta::invmat: ... must be a symmetric matrix with a positive definite imaginary part – ssch Feb 14 '13 at 12:45
• The function is defined for all $\Omega$ with positive definite imaginary part, just as the error message says. In my case I'm only interested in the region $(x,y,z)\in[0,1]\times [0,1-x]\times [0,1]$. Also I'm only interested in the case where the second argument is $s=(0,0)$, and the reason I introduced a non-zero value for $s$ is only so it will allow me to take the derivative. – Emil Bostrom Feb 14 '13 at 14:20
• Wow, this was a really tricky function to do anything useful with, seems it evaluates really slowly at some points, doing an interpolation of it now I'll give it another 15minutes :p – ssch Feb 14 '13 at 16:53

You could try a numeric approximation for the second derivatives along the lines of:

o[x_, y_, z_] := {{(1/z - 1)^(1/2)*I/Pi*(1 - y),      (1/z - 1)^(1/2)* I/Pi*(-1 + x + y)},
{(1/z - 1)^(1/2)*I/Pi*(-1 + x + y), (1/z - 1)^(1/2)*I/Pi*(1 - x)}};
st = SiegelTheta;
f2 = Compile[{{x, _Complex}, {y, _Complex}, {z, _Complex}, {h, _Complex}},
(-st[o[x, y, z], {0, 2 h}] + 16 st[o[x, y, z], {0, h}] -
30 st[o[x, y, z], {0, 0}] + 16 st[o[x, y, z], {0, -h}] -
st[o[x, y, z], {0, -2 h}])/(12 h h)];


The problem is that NIntegrate has difficulties dealing with

NIntegrate[Chop[N[f2[x, y, z, 1/100.]]], {x, 0, 1}, {y, 0, 1 - x}, {z, 0, 1}]


Because the function is ill behaved when z is near 0 around {0,0}. Look:

ListPlot3D[Flatten[Chop@ Table[{x, y, f1[x, y, .01, 1/100]}, {x, 0.01, .99, .1},
{y, 0.01, 1 - x, .1}], 1], PlotRange -> Full] so you'll need to find a suitable integration method (I haven't tried hard, but couldn't find it)

• I like the approximative approach and will definitely try that out. Ideally however I would like to have the function in an exact (or close to) form since I'm going to use it in quite precise calculations. I've tried to compile the function starting from the definition according to Compile[{x, y, z}, Sum[n^2*Exp[-n^2*(1/z - 1)^\[Alpha]*(1 - y) - m^2 (1/z - 1)^\[Alpha]*(1 - x) - 2*n*m (1/z - 1)^\[Alpha]*(-1 + x + y)], {n, -Infinity, Infinity}, {m, -Infinity, Infinity}]] but I'm having trouble integrating that as well. Any thoughts? – Emil Bostrom Feb 14 '13 at 19:08
• By the way Alpha above should be 1/2. – Emil Bostrom Feb 14 '13 at 19:19
• @EmilBostrom already tried that to no avail – Dr. belisarius Feb 14 '13 at 20:36
• I can get it to work if I compile using NSum instead of Sum, so the problem must have to do with Mathematica trying to do the sum symbolically. – Emil Bostrom Feb 15 '13 at 11:02