# Tag Info

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

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

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 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 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 ... 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 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 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 You could triangulate the domain (with say a Delaunay method), and then integrate over the triangles using standard integration techniques. This would be a finite element method, but just for integrating a function rather than solving a differential equation. 5 To give you a better picture, what I mean with my comment, consider Stokes Theorem$$\int_D db^{n-1}=\int_{\partial D}b^{n-1}$$be b^{n-1} an arbitrary 1-form$$b^{1}= a_1 dx + a_2 dy$$leads to$$\int_D f(x,y)dxdy = \int_D\left(\frac{\partial a_2}{\partial x}-\frac{\partial a_1}{\partial y}\right)dxdy$$Without loss of generality, a_1=0 one needs to ... 5 Replace the last line by plot(X, [F(x)[0] for x in X]) That should do it. Edit: you can define your function F as def F(x): try: return [integrate.quad(f, 0, y) for y in x] except TypeError: return integrate.quad(f, 0, x) 5 If you know the Chebyshev expansion for f(z), why don't you formally integrate the polynomials using the recurrence relation for Chebyshev polynomials ? The Clenshaw-Curtis method is based on this approach (combined with an intelligent use of the FFT). 5 With the substitution y=x+t, the t-dependency can be factored out:$$I(t)=\int_{-\infty}^\infty e^{-(y-t)}F(y)\;dy = e^t \int_{-\infty}^\infty e^{-y}F(y)\;dy = C\cdot e^t This still leaves you with the problem of computing $C$, but it is merely a constant number and the $t$-dependency is entirely described by the factor $e^t$. The two-sided Laplace ...

5

Step 1: analyze the recurrence I implemented the recurrence in Python, using basic numpy and scipy.special for the erf function. Why? Because it is simple and cheap (for reasonable $p$). import numpy as np import scipy.special as sp def F1(z): return np.sqrt(np.pi/2)*(1.0+sp.erf(z/np.sqrt(2))) def F2(z): return F1(z)*z + np.exp(-z**2/2.0) def F(...

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