4
$\begingroup$

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.

Improve this question
$\endgroup$
6
  • $\begingroup$ Can you post the Mathematica expression for Ω $\endgroup$
    – ssch
    Commented Feb 14, 2013 at 10:22
  • $\begingroup$ 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)}} $\endgroup$
    – Emil Bostrom
    Commented Feb 14, 2013 at 10:52
  • $\begingroup$ 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 $\endgroup$
    – ssch
    Commented Feb 14, 2013 at 12:45
  • $\begingroup$ 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. $\endgroup$
    – Emil Bostrom
    Commented Feb 14, 2013 at 14:20
  • $\begingroup$ 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 $\endgroup$
    – ssch
    Commented Feb 14, 2013 at 16:53

1 Answer 1

Reset to default
3
$\begingroup$

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)

Improve this answer
$\endgroup$
4
  • $\begingroup$ 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? $\endgroup$
    – Emil Bostrom
    Commented Feb 14, 2013 at 19:08
  • $\begingroup$ By the way Alpha above should be 1/2. $\endgroup$
    – Emil Bostrom
    Commented Feb 14, 2013 at 19:19
  • $\begingroup$ @EmilBostrom already tried that to no avail $\endgroup$
    – Dr. belisarius
    Commented Feb 14, 2013 at 20:36
  • $\begingroup$ 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. $\endgroup$
    – Emil Bostrom
    Commented Feb 15, 2013 at 11:02

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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