# Numerical integration with singularities

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.

• How big are $j,N,A(m,n)$, what are the typical values? Can you post a small reproducible example of parameter values for which the problem occurs? How is the integral singular, the tan terms don't look singular to me? A likely source of failure is the fast oscillation due to $e^{iy(n-m)}$, which adaptive quadrature routines like dblquad would usually struggle with. – Kirill Sep 10 '15 at 14:28
• I've made some aclarations to the text. The integral itself is not singular, but the quantities involved in the calculus are very big and produce overflow in the calculus. Maybe simply solving the oscillation the dblquad can make something useful with that. – VictorSeven Sep 10 '15 at 14:53
• $x^{2j+m+n}/(1+x^2)^{2j}$ as a function of $x=\tan\frac{x}{2}\geq 0$ is always between $0$ and $1$. When $x$ is large, you can evaluate it as $x^{m+n-2j}/(x^{-2}+1)^{2j}$, which should address the overflow issues. – Kirill Sep 10 '15 at 14:56

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

• Brilliant. That idea worked perfect for the overflows. Also dividing the double integral using two different methods -and the idea that the integral in $y$ would need only a few iterations- works perfect. I did the integral with j=60, 20 times, only in 90 s with C++. Thank you! :D – VictorSeven Sep 11 '15 at 10:31
• @V_Programmer Glad it worked. – Kirill Sep 11 '15 at 13:17