6
$\begingroup$

I would like to perform the numerical integration of an integral of the form $$ \int_{-\infty}^\infty e^{i \omega 0+} G(i \omega, \mathbf{v}) d \omega ,$$ or, using the symmetry $G(i\omega)^* = G(-i \omega)$, $$ \int_0^\infty \Re(e^{i \omega 0+} G(i \omega, \mathbf{v}) ) d\omega ,$$ where $\mathbf{v}$ is a vector of additional parameters. $G(i \omega)$ is a rational function and behaves for large $\omega$ as $\omega^5/\omega^6 = 1/\omega$ and it can not be integrated analytically.

I split the integral, i.e. $\int_0^\infty = \int_0^a + \int_a^\infty$, where $a$ is a large positive constant, such that the approximation $G(i \omega) \approx 1/i \omega$ holds for $\omega \geq a$. The infinite integral gives an additive constant. The convergence factor in $\int_0^a$ can now be dropped because it is a finite integral.

Now there is this problem: I use the quad() function from SciPy to evaluate the finite integral for set of different $\mathbf{v}$. If I set the upper limit in quad() as $a=10^8, 10^9$ or $a = \infty$ then it gives consistent and correct results. But for $a = 10^{10} - 10^{14}$, the results are rubbish. For $a = \infty$ quad() gives a warning message that the integrand may be slowly convergent and that there is a roundoff error detected in the extrapolation table.

So obviously there is a numerical problem and I don't know what quad() exactly does, when one enters $\infty$ as a limit. From the source code it seems that quad() maps the big finite interval onto $[0,1]$, rescales the integrand and uses a Gauss-Kronrod method for integration, but it is not clear (to me) how big this scaling respectively the finite interval is.

Are there other reliable methods to perform a numerical integration over large intervals? Is a method where one divides the big interval into a number of subintervals with fixed width and subsequent integration (e.g. with Simpson's rule) recommended?

$\endgroup$
3
$\begingroup$

Ok, I fiddled around with it and think I got a workaround. In fig. 1 we see a graph of $G(i \omega)$ for a typical set of parameters (printed on the graph, you can ignore them) and realize that $G(i \omega)$ falls off as $1/\omega^2$ for $\omega \geq 10^6$. So the biggest contribution to the integral comes from the interval $[0, 10^6]$. Here is a list of values for the integral that I get for this function when I use $\tt{scipy.integrate.quad()}$ with different values for the upper limit $a$:

  • $a$, result
  • $10^8$, $1.54368741837$
  • $10^9$, $1.54428315319$
  • $10^{10}$, $1.54434272669$
  • $10^{11}$, $-6.61927531679 \cdot 10^{-7}$
  • $10^{12}$, $-6.61927690454 \cdot 10^{-8}$
  • $10^{14}$, $-6.61927691368 \cdot 10^{-10}$
  • numpy.inf, $1.54435019948$

I suppose that the rubbish results starting at $10^{11}$ arise because the $10^6$ is only a very very small percentage of $10^{11}$ and the quad()-algorithm does not sample enough values of $G$ in the essential region of the interval and therefore gives a severly underestimated value of the integral. This problem is somehow circumvented when one uses numpy.inf as an upper limit. But I want to avoid the use of it, because I do not really understand how quad() handles this case.

The solution to this problem is the following. We rewrite the integral as an initial value problem (IVP) of a first order ODE: \begin{align*} I(x) &:= \int_0^x~G(i\omega, \mathbf{v}) d\omega \\ \Rightarrow \frac{dI}{dx} &= G(i x, \mathbf{v}), \; I(0) = 0, \end{align*} and we are interested in the value $I(a)$ for $a$ very large. Now we can use an ODE integrator to solve this IVP. I used a Runge-Kutta Cash-Karp method (with adaptive stepsize!) that I translated from the book "Numerical Recipes" (2nd ed.) and the standard ODE solver $\tt{scipy.integrate.odeint()}$. Here are the results for $10^8 - 10^{14}$ in this order.

Runge-Kutta Cash-Karp (RKCK): 1.54368746, 1.5442832 , 1.54434278, 1.54434874, 1.54434934, 1.5443494

odeint: 1.54368718, 1.54428289, 1.54434245, 1.5443484 , 1.54434901, 1.54434908

RKCK needs around 100 steps on average to solve the integral, odeint() more than 200 steps. Using an ODE solver seems to resolve the problem because it can adapt its stepsize and sample enough relevant function values to give the correct result of the integral.

The graph of $G(i \omega)$ for a typical set of parameters

$\endgroup$
3
$\begingroup$

SciPy's integration tools are based on the QUADPACK numerical integration library. If you really want to get into the gory details of how quad() treats improper integrals and what specific error messages mean, you can check out the source code and documentation at:

http://www.netlib.org/quadpack/

The specific routine that is called for infinite intervals is QAGI, and helpful documentation is available in the doc file. If Fortran isn't your thing, you could also look at the equivalent C functions in the GSL:

https://www.gnu.org/software/gsl/manual/html_node/Numerical-Integration.html

As for your specific problem, I think your interpretation is correct. The adaptive quadrature works by sampling the function at 22 points, fitting a polynomial, and integrating that. This result is compared to the lower-order polynomial resulting from only 15 of those points. If the two are in close agreement, then integration is halted. Otherwise, the domain is automatically bisected and the pieces evaluated separately.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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