20

In one dimension, you can map your infinite interval to a finite interval using integration by substitution, e.g. $$\int_a^b f(x)\,\mathrm dx \quad=\quad \int_{u^{-1}(a)}^{u^{-1}(b)}f(u(t))u'(t)\,\mathrm dt$$ Where $u(x)$ is some function that goes off to infinity in some finite range, e.g. $\tan(x)$: $$\int_{-\infty}^{\infty}f(x)\,\mathrm dx \quad=\quad ...


20

I think this is not quite what you had in mind, but for the sake of completeness, let's start with some basics. Most quadrature formulas such as Newton-Cotes and Gauss are based on the idea that in order to evaluate the integral of a function approximately, you can approximate the function by, e.g., a polynomial that you can then integrate exactly: $$ \int_a^...


19

I'm not entirely familiar with what's now done for cubatures (multidimensional integration), so I'll restrict myself to quadrature formulae. There are a number of effective methods for the quadrature of oscillatory integrals. There are methods suited for finite oscillatory integrals, and there are methods for infinite oscillatory integrals. For infinite ...


15

The question has two very different subquestions. I will address the first one only. Matlab's version runs on average 24 times faster than my python equivalent! The second one is subjective. I would say that letting know the user that there is some problem with the integral is a good thing and this SciPy behavior outperforms the Matlab`s one to keep it ...


12

I've written my own integrator, quadcc, which copes substantially better than the Matlab integrators with singularities, and provides a more reliable error estimate. To use it for your problem, I did the following: >> lambda = 0.00313; kappa = 0.00825; nu = 0.33; >> x = 10; >> E = @(r) r.^4.*(lambda*sqrt(kappa^2 + r.^2)).^(-nu-5/2) .* ...


12

Use Plancherel's theorem to evaluate this integral. The basic idea is that for two functions $f,g$, $$ I=\int_{-\infty}^{\infty} f(x) g^*(x)dx = \int_{-\infty}^{\infty} F(k) G^*(k) dk $$ where $F,G$ are the Fourier transforms of $f,g$. Your functions both have relatively small support in the spectral domain. Here, $\sin x / x \rightarrow \text{rect}(k)$ ...


11

The standard way of doing it is to extract from the expression for $f(x)$ an exponential prefactor, transform that to $e^{-x^2}$, and then use Gaussian quadrature rules (or Gauss Kronrod) with this as a weight. If $f$ is smooth, this usually gives excellent results. In $R^3$, the same works with weight $e^{-|x|^2}$, and appropriate cubature formulas can be ...


11

Besides "multidimensional vs. single-dimensional" and "finite range vs. infinite range", an important categorization for methods is "one specific type of oscillator (usually Fourier-type: $\sin(t)$, $\exp(it)$, etc, or Bessel-type: $J_0(t)$, etc.) vs. more general oscillator ($\exp(i g(t))$ or even more general oscillators $w(t)$)". At first, oscillatory ...


11

As Pedro points out, Levin-type methods are the best established methods for these kinds of problems. Do you have access to Mathematica? For this problem, Mathematica will detect and use them by default: In[1]:= e[r_] := r^4 (l Sqrt[k^2 + r^2])^(-v - 5/2) BesselK[-v - 5/2, l Sqrt[k^2 + r^2]] In[2]:= {l, k, v} = {0.00313, 0.00825, 0.33}; In[3]:= Block[{...


11

First of all, you need to ask yourself the question if you need an all-round quadrature routine that should take an integrand as a black box. If so, you cannot but go for adaptive quadrature where you hope that the adaptivity will catch "difficult" spots in the integrand. And that is one of the reasons Piessens et al. chose for a Gauss-Kronrod rule (this ...


10

I am the author of CQUAD in the GSL. The interface is almost identical to that of QAGS, so if you've used the latter, it should not be difficult at all to try the former. Just remember not to convert your NaNs and Infs to zeros in the integrand -- the code will deal with these itself. The routine is also available in Octave as quadcc, and in Matlab here. ...


10

This can be done by integration by parts: $$ \int^\infty_1 x e^{-ax} = \frac{-1}{a} x e^{-ax}\mid^\infty_1 - \frac{-1}{a} \int^\infty_1 e^{-ax} = \frac{e^{-a}}{a} + \frac{e^{-a}}{a^2} = \frac{a+1}{a^2} e^{-a} $$ and continuing on by induction $$ \int^\infty_1 x^k e^{-ax} = \frac{-1}{a} x^k e^{-ax}\mid^\infty_1 - \frac{-k}{a} \int^\infty_1 x^{k-1} e^{-...


10

There are a bunch of pitfalls when it comes to tanh-sinh quadrature, one being that the integrand needs to be evaluated very closely to the interval boundaries, at distances less than machine precision, e.g., 1.0 - 1.0e-20 in the original example. When this point is evaluated, it rounds to 1.0 at which f has a singularity, and anything can happen. You're ...


10

You can evaluate this using the Björck-Pereyra algorithm for solving Vandermonde systems, because you are evaluating $b^\top V^{-1}$ with $b=(2,0,\frac23,0,\frac25,0,\ldots)$, and the algorithm is known to be forward-stable (see Error analysis of the Björck-Pereyra algorithms for solving Vandermonde systems by Nick Higham, http://www.maths.manchester.ac.uk/~...


9

This is a linear integral equation (because the right hand side is linear in $h$) and presumably the problem is ill-posed (for one, because of the multiplication by $\sin(\theta)$ but also because the "integral kernel" is smooth). Methods to solve this (approximately) are general Galerkin methods or collocation methods. For Galerkin methods you take a ...


9

For large $x$, you can use the same approach as the one in the paper you cite. They define these auxiliary functions that link $\mathrm{Si},\mathrm{Ci}$ with $\sin,\cos$, and the point of why this works is that the relationship between these two pairs of functions is in general a $2\times 2$ matrix, but this matrix only has two independent functions in its ...


8

You should be able to get accurate results with mpmath, a Python module for arbitrary-precision floating-point computations. There are examples of integration with singularities in the documentation. You'll want to explicitly tell it to break up the interval: from mpmath import * f = lambda x,y,z: 1./(x**2+y**2+z**2)**1./3 quad(f,[-1,0,1],[-1,0,1],[-1,0,1]...


8

You're integrating a function on a 2-dimensional manifold embedded in $\mathbb{R}^{3}$; books in analysis on manifolds (like Munkres' accessible book, or Lee's books on manifolds) are helpful in discussing the theory defining this type of integral. Let's suppose that $f$ is a real-valued function defined on the manifold $\mathcal{M}$, which is your 4-node 3-...


8

There are standard methods for these types of quadrature in Python, in NumPy and SciPy: Gauss-Laguerre quadrature Gauss-Legendre quadrature Gauss-Hermite quadrature (as noted in your post) Gauss-Chebyshev quadrature QUADPACK adaptive quadrature adaptive Gaussian quadrature and other routines. Given a degree of quadrature (essentially, the degree of the ...


8

You may use the ideas of error extrapolation as one uses it to construct high-order Runge Kutta methods. Depending on the function that you interpolate, the interpolation error $I - I_h$, where $I$ is the actual integral value and $I_h$ the value obtained by the piecewise trapezoidal rule, may have a smooth asymptotic expansion $$ I - I_h = C_p h^p + C_{p+...


7

Disclaimer: I wrote my PhD thesis on adaptive quadrature, so this answer will be severely biased towards my own work. GSL's QAGS is the old QUADPACK integrator, and it is not entirely robust, especially in the presence of singularities. This usually leads to users requesting far more digits of accuracy than they actually need, thus making the integration ...


7

Is it important to have the most efficient rules possible? Do you require all weights to be positive and/or all points to be inside the region? Two general resources: The Encyclopaedia of Cubature Formulas and referenced papers. The Libmesh source contains good choices of cubature rules for all orders, with and without negative weights, thanks to John ...


7

For one-dimensional quadrature, you can check the book on Quadpack (a golden oldie but still very relevant in one-dimensional quadrature) and the techniques used in the algorithm QAGI, an automatic integrator for an infinite range. Another technique is the double-exponential quadrature formula, nicely implemented for an infinite interval by Ooura. For ...


7

It depends on how big $m$ and $n$ are (assuming that $x'$ ranges from $0$ to $a$, though the notation suggests the opposite if read literally). You should probably split the intervals into at least $m$ resp $n$ subintervals, integrate the pieces separately, and sum the results. If $f$ is smooth, a Gaussian rule is probably more accurate for the same effort. ...


7

Using product quadrature rules for multi-dimensional integrals suffers from the so-called curse of dimensionality. An $O(N^{-2})$ accurate rule using N evaluations in one-dimension is generally $O(N^{-2})$ accurate when applied as a product rule in multiple dimensions, but there will be $M = N^d$ evaluations required. So the accuracy is $O(M^{-2/d}).$ The ...


7

The problem is almost definitely with how QUADPACK (which is the backend used by scipy.integrate.quad) handles numerically small integrands. Essentially the integrand is so small (at $x=0$ it is $6.58\times 10^{-12}$), that it is comparable to the absolute error tolerance. I'm not sure quite why it says it's "probably divergent", but that doesn't matter. I ...


7

The key to the evaluation of oscillatory integrals is to truncate integral at the right point. For this example you need to choose upper limit of the form $$ \pi\mathbb{N}+\frac{\pi}{2} $$ Before explaining why it should work, let me at first show that it actually produces good results. Asymptotics It is easy to guess that asymptotic series has the ...


6

If you remember how quadrature worked, then you'll also know one way how to approximate infinite integrals. Namely: for quadrature, you approximate the function $f(x)$ you want to integrate by something similar, say a polynomial $\tilde f(x)$ (or a piecewise polynomial) for which you can write down the integral analytically. You get $\tilde f$ from $f$ by ...


6

There are a number of "corrected" integration rules which invoke derivatives of the endpoints. One simple example is the corrected trapezoidal rule. Suppose we wish to approximate the integral $$ \int_a^bf(x)\,dx. $$ Let $n$ be an integer, and $h=(b-a)/n$. Then the trapezoidal rule $$ T = \frac{h}{2}\bigl(f(a)+2f(a+h)+2f(a+2h)+\cdots+2f(a+(n-1)h)+f(b)\...


6

Since the outer integral $\int_0^{2\pi} F(\theta) d\theta$ has a periodic integrand, a trapezoidal rule should work nicely, picking up extra accuracy when $F(\theta)$ is smooth. For the inner integral you want to take advantage of the known "weight function" $r$ that appears in the integrand, $F(\theta) = \int_0^{r(\theta)} f(r,\theta) r\; dr$. One way to ...


Only top voted, non community-wiki answers of a minimum length are eligible