Numerical integration with singularities

Question

I need to compute some integrals numerically. The integrand is this:

$$f(x,y) = \left ( \sum_{mn=-j}^{j}A(m,n)\dfrac{\tan^{2j+m+n}(x/2)}{(1+\tan^2(x/2))^{2j}}e^{iy(n-m)} \right )^{N}$$

Note: sums are finite. And I have to compute

$$\int_0^\pi \int_0^{2\pi}f(x,y)dydx$$

I have to do this for different values of $N$. I was working in C++ with a simple integrator, but the problem is that the limits in $x$ are $[0,\pi]$. When $j$ is big and $x\rightarrow\pi$, I have overflows/underflows everywhere.

I've also tried to use Python/Scipy, with the dblquad function. However, it takes forever to compute the integral. The results aren't good. I have to compute it more or less a hundred of times, so I need something reasonably fast.

I'm really stuck with this. Any idea to get a numerical result from the integral is well received.

EDIT: Common numerical values as requested.

$$j\in[20,150]$$ $$n,m\in[-j,j]$$ $$N=1,2,...,5$$

That means that we have up to ~$\tan^{4jN}(x/2)$ both in numerator/denominator. When $x$~$\pi$ that becomes huge! I've obtained a lot of NaNs in C++. $A(m,n)$ is simply a coefficient obtained from some non-zero constants.

Important note: althought the term $e^{iy(n-m)}$ is there, $f(x,y)\in\mathbb{R}$. Imaginary parts always dissapears.

Answer 1

I think there are two things that could be done here.

First, the way you evaluate $$ s(x) = \frac{\tan^{2j+m+n}\frac x2}{(1+\tan^2\frac{x}{2})^{2j}} $$ seems to be unstable, even though we always have $0\leq s(x)\leq 1$. Thus the overflows can be addressed by rewriting this expression as $$ \frac{z^{2j+m+n}}{(1+z^2)^{2j}} = \frac{(1/z)^{2j-m-n}}{(1+z^{-2})^{2j}}, $$ and using this alternative expression when $z = \tan\frac{x}{2} > 1$.

Second, the oscillatory nature of the integral will lead to trouble when trying to integrate it using an adaptive quadrature routine like scipy's quad (which uses QUADPACK). The problem is that the integration error for the integrand $f(y)$ will behave as some high derivative $f^{(m)}$, which for an oscillatory function $e^{i y p}$ will be large, proportional to $p^m$.

I think the way to address this is to use trapezoid rule, which is very accurate when applied to periodic functions. I think an appropriate method might be to write the integral as $$ \approx \int_0^\pi \left( 2\pi N_y^{-1}\sum_{k=0}^{N_y-1} f(x, 2\pi k/N_y) \right)\,dx, $$ and evaluate the remaining one-dimensional integral using quad; a suitable value of $N_y$ should be found experimentally/adaptively, it should be quite small - I'm only writing it this way because it seems scipy doesn't already include an adaptive trapezoid quadrature function. The double-exponential rule (which isn't used by quad) is also a good candidate to try here for the outer integral.

Finally, the third thing to try is the change of variable $t = \tan\frac{x}{2}$, with $$ \frac{2}{1+t^2}\,dt = dx, \qquad \int_0^\pi (\cdot)dx = \int_0^\infty (\cdot)\frac{2 dt}{1+t^2}. $$ This might make it easier for quad to deal with the function, though I can't be sure without checking.

Here's the code I wrote in Mathematica to make sure any of this works:

Clear[a, f, x2, g1, g2]
a[m_, n_] := 1/(m^2 + n^2 + 1)
x2[j_, s_, t_?NumericQ] := 
 If[t > 1, (1/t)^(2 j - s)/(1 + (1/t)^2)^(2 j), t^(2 j + s)/(1 + t^2)^(2 j)]
f[j_, q_, x_, y_] := 
 Sum[a[m, n] x2[j, m + n, Tan[x/2]] Exp[I y (m - n)], {m, -j, j}, {n, -j, j}]^q // Re
g1[j_, q_] := NIntegrate[f[j, q, x, y], {x, 0, \[Pi]}, {y, 0, 2 \[Pi]}]
g2[j_, q_] := Module[{g},
  g[x_?NumericQ] := 
   NIntegrate[f[j, q, x, y], {y, 0, 2 \[Pi]}, Method -> "TrapezoidalRule"];
  NIntegrate[g[x], {x, 0, \[Pi]}, Method -> "DoubleExponential"]
  ]

It evaluated $j=20, N=2$ as 3.988585896935342e-6 in about 26s, with relative error $7.9\times 10^{-8}$ (checked with extra precision).