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


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

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

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


11

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

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

First of all, your function $x\sin(\frac{1}{x})$ is singular in $x=0$. You might want to add an if clause like this: def f(x): if abs(x) < 1e-10: res = x else: res = x*sin(1/x) but this does hurt speed (masked arrays would be better). The reason why your code doesn't work is because scipy.quad only accepts a single value as ...


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

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

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

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


6

The issue here is that the error term increases very quickly as $s$ increases. The error term is, according to this article, is bounded by $$ |E| < \frac{n!^2}{(2n)!}\max |f^{(2n)}(t)|. $$ Now, asymptotically $$ \frac{n!^2}{(2n)!} \approx \frac{\sqrt{\pi n}}{4^n}, $$ and $$ |f^{(2n)}(t)| = s^3(1-s^2)^{2n} e^{-s^2t} < s^3(1-s^2)^{2n}, $$ meaning that ...


6

Yes. As you may know, numerical differentiation and integration is closely related to (polynomial) interpolation: The idea to approximately differentiate or integrate a given function is to approximate it with a function (often an interpolating polynomial) that can be differentiated or integrated exactly. For example, the standard central difference quotient ...


6

The point of @WolfgangBangerth is exactly what I mentionend in my comment, so I'd always try this first. In the best case, with millions of partitions $[a_{i},a_{i+1}]_{i\in \{0,\ldots,N-1\}}$ (where $a_{i} < a_{i+1}$), the length of these is quite small such that some low-order method could be appropriate. This will require roughly $N \times (s-1)$ ...


6

Ooura's method for Fourier sine integrals works here, see: Ooura, Takuya, and Masatake Mori, A robust double exponential formula for Fourier-type integrals. Journal of computational and applied mathematics 112.1-2 (1999): 229-241. I wrote an implementation of this algorithm but never put in the work to get it fast (by, say caching nodes/weights), but ...


5

Physically, this problem arises when the two input nodes are "impossibly far apart" according to the planar wavefront approximation. Specifically, if two nodes separated by distance $\Delta x \sqrt{2}$ on the grid have $T$ values differing by more than $\Delta x \sqrt{2}$, there is no solution plane that fits both distance values. Even if the quadratic ...


5

Apart from the methods based on Newton-Cotes mentioned in the other answers, there is what is now called Gauss-Turán quadrature (see e.g this and this, as well as the venerable reference by Walter Gautschi). This is a generalization of the usual Gaussian quadrature, where one can now exploit the knowledge of a function's derivatives in finding an optimal set ...


5

If you don't have access to Mathematica, you could write a Levin-type (or other specialized oscillatory) method in Matlab as Pedro suggests. Do you use the chebfun library for Matlab? I just learned it contains an implementation of a basic Levin-type method here. The implementation is written by Olver (one of the experts in the oscillatory quadrature field)....


5

Your approach to selecting $h$ based on $N$ isn't a very sophisticated one and probably isn't getting $h$ to be small enough. See these notes for a more reasonable way to adjust $h$ and estimate the error in the integral evaluation: http://www.davidhbailey.com/dhbpapers/dhb-tanh-sinh.pdf


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