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


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


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

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


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

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[{...


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

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'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

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


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

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

First of all, from the first paragraph of your attempts at a solution, I assume that the $z_j$ are non-negative? In that case, the integrand has no real problematic points (it's monotonous, decreasing). The only difficulty is the infinite integration range. And then it all depends on the accuracy you want... Do you want machine epsilon precision (or there ...


7

Instead of directly integrating over the area, it is often more convenient to use the divergence theorem to replace the area integral with an integral over the boundary edges. The divergence theorem in three dimensions for a vector function $F$ can be written $$ \int_V \nabla \cdot {\bf F} dv = \int_S {\bf F} \cdot ds $$ That is, an integral over the ...


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

You could just go for the direct final answer of the integral: $$\int_{a}^{b}\vert {\sin(x)}\vert dx = 2\left(\left\lfloor\frac{b}{\pi}\right\rfloor - \left\lfloor\frac{a}{\pi}\right\rfloor\right) + \cos\left(\left\{\frac{a}{\pi}\right\}\pi\right) - \cos\left(\left\{\frac{b}{\pi}\right\}\pi\right)$$ where $\left\lfloor c\right\rfloor$ denotes the flooring ...


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

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

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


6

You can't bound the accuracy of any numerical integration scheme without some control over how crazy the integrand can be. Smoothness is obviously helpful, but if the derivatives aren't bounded, then it's possible to introduce a $C^{\infty}$ bump into the integrand that can have an arbitrarily large impact on the value of the integral and you might not ever ...


6

For the Hankel transform, one can classify the methods into four major groups: Numerical quadrature-based. Fourier-based ones. Asymptotic expansion of Bessel into sines and cosines. Projection-slice methods. The following paper gives a nice overview of these methods (types of methods): M. J. Cree and P. J. Bones, "Algorithms to numerically evaluate the ...


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