Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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.

share|improve this question
    
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
show 1 more comment

1 Answer

up vote 0 down vote accepted

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]

Mathematica graphics

so you'll need to find a suitable integration method (I haven't tried hard, but couldn't find it)

share|improve this answer
    
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 –  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
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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