When the interval of integration or the integrand itself is unbounded, we say an integral is improper . Improper integrals present particular challenges to numerical computation.
9.7.1 Infinite interval ¶ When the integration domain is ( − ∞ , ∞ ) (-\infty,\infty) ( − ∞ , ∞ ) , the integrand has to decay as x → ± ∞ x \to \pm \infty x → ± ∞ in order for the improper integral to be finite. This fact brings up the possibility of truncating the domain:
∫ − ∞ ∞ f ( x ) d x ≈ ∫ − M M f ( x ) d x . \int_{-\infty}^{\infty} f(x)\, dx \approx \int_{-M}^{M} f(x)\, dx. ∫ − ∞ ∞ f ( x ) d x ≈ ∫ − M M f ( x ) d x . This integral can be discretized finitely by, say, the trapezoid formula, or an adaptive integrator. However, this approach can be inefficient.
Consider the integral of f ( x ) = 1 / ( 1 + x 2 ) f(x)=1/(1+x^2) f ( x ) = 1/ ( 1 + x 2 ) ,
∫ − ∞ ∞ 1 1 + x 2 d x = π . \int_{-\infty}^\infty \frac{1}{1+x^{2}}\, dx = \pi. ∫ − ∞ ∞ 1 + x 2 1 d x = π . In this case we can easily estimate the effect of truncation on the result. For large M M M ,
∫ M ∞ f d x ≈ ∫ M ∞ x − 2 d x = M − 1 . \int_M^\infty f\,dx \approx \int_M^\infty x^{-2}\,dx = M^{-1}. ∫ M ∞ f d x ≈ ∫ M ∞ x − 2 d x = M − 1 . The same estimate applies to the integral over ( − ∞ , − M ) (-\infty,-M) ( − ∞ , − M ) . To get eight digits of accuracy, for instance, we need to truncate with M > 2 × 1 0 8 M > 2 \times 10^8 M > 2 × 1 0 8 .
In order to do better than direct truncation, we want to encourage the function to decay faster. In practice this means a change of variable, x ( t ) x(t) x ( t ) . If ∣ x ( t ) ∣ |x(t)| ∣ x ( t ) ∣ grows rapidly as ∣ t ∣ → ∞ |t| \to \infty ∣ t ∣ → ∞ , then f ( x ( t ) ) f(x(t)) f ( x ( t )) will decay more rapidly in t t t than in x x x .
One way to accomplish this feat is to use
x ( t ) = sinh ( sinh t ) . x(t) = \sinh\left( \sinh t \right). x ( t ) = sinh ( sinh t ) . Noting the asymptotic behavior as t → ± ∞ t \rightarrow \pm\infty t → ± ∞ that
∣ sinh ( t ) ∣ ∼ 1 2 e ∣ t ∣ , \left| \sinh(t) \right| \sim \frac{1}{2} e^{ |t| }, ∣ sinh ( t ) ∣ ∼ 2 1 e ∣ t ∣ , we find that in the same limits,
x ( t ) ≈ ± 1 2 exp ( 1 2 e ∣ t ∣ ) . x(t) \approx \pm \frac{1}{2} \exp\left( \frac{1}{2} e^{ |t| } \right). x ( t ) ≈ ± 2 1 exp ( 2 1 e ∣ t ∣ ) . Thus, (9.7.4) is often referred to as a double exponential transformation.
By the chain rule,
∫ − ∞ ∞ f ( x ) d x = ∫ − ∞ ∞ f ( x ( t ) ) d x d t d t = ∫ − ∞ ∞ f ( x ( t ) ) cosh ( sinh t ) cosh t d t . \begin{split}
\int_{-\infty}^\infty f(x)\, dx &= \int_{-\infty}^\infty f(x(t))\frac{dx}{dt}\, dt \\
&= \int_{-\infty}^\infty f(x(t))\, \cosh\left( \sinh t \right) \cosh t \, dt.
\end{split} ∫ − ∞ ∞ f ( x ) d x = ∫ − ∞ ∞ f ( x ( t )) d t d x d t = ∫ − ∞ ∞ f ( x ( t )) cosh ( sinh t ) cosh t d t . The exponential terms introduced by the chain rule grow double exponentially, but the more rapid decay of f f f in the new variable more than makes up for this.
Example 9.7.2 (Decay by transformation)
Consider again f ( x ) = 1 / ( 1 + x 2 ) f(x)=1/(1+x^2) f ( x ) = 1/ ( 1 + x 2 ) from Example 9.7.1 , with x ( t ) x(t) x ( t ) given by (9.7.4) . As t → ∞ t\to\infty t → ∞ ,
f ( x ( t ) ) ≈ x − 2 ≈ 4 exp ( − e t ) . f(x(t)) \approx x^{-2} \approx 4 \exp\left( -e^{t} \right). f ( x ( t )) ≈ x − 2 ≈ 4 exp ( − e t ) . The chain rule terms in (9.7.7) become
cosh ( sinh t ) cosh t ≈ 1 2 exp ( 1 2 e t ) ⋅ 1 2 e t , \cosh\left( \sinh t \right) \cosh t
\approx \frac{1}{2} \exp\left( \frac{1}{2} e^t \right) \cdot \frac{1}{2} e^t, cosh ( sinh t ) cosh t ≈ 2 1 exp ( 2 1 e t ) ⋅ 2 1 e t , yielding a product that is roughly
2 exp ( − 1 2 e t ) . 2 \exp\left( -\frac{1}{2} e^t \right). 2 exp ( − 2 1 e t ) . The total integrand in (9.7.7) therefore has double exponential decay in t t t , essentially because of the squaring of x x x in the denominator of f f f . The same result holds as t → − ∞ t\to-\infty t → − ∞ .
Example 9.7.2 using Plots
f(x) = 1 / (1 + x^2)
plot(f, -4, 4, layout=(2, 1),
xlabel=L"x",
yaxis=(:log10, L"f(x)", (1e-16, 2)),
title="Original integrand")
ξ(t) = sinh( π * sinh(t) / 2 )
dξ_dt(t) = π/2 * cosh(t) * cosh(π * sinh(t) / 2)
g(t) = f(ξ(t)) * dξ_dt(t)
plot!(g,-4, 4, subplot=2,
xlabel=L"t",
yaxis=(:log10, L"f(x(t))\cdot x'(t)", (1e-16, 2)),
title="Transformed integrand")
This graph suggests that we capture all of the integrand values that are larger than machine epsilon by integrating in t t t from -4 to 4.
Example 9.7.2 f = @(x) 1 ./ (1 + x.^2);
clf, subplot(2, 1, 1)
fplot(f, [-4, 4]); set(gca, 'yscale', 'log')
xlabel('x'), ylabel('f(x)'), ylim([1e-20, 1])
title('Original integrand')
x = @(t) sinh( pi * sinh(t) / 2 );
chain = @(t) pi/2 * cosh(t) .* cosh( pi * sinh(t) / 2 );
integrand = @(t) f(x(t)) .* chain(t);
subplot(2, 1, 2)
fplot(integrand, [-4, 4]); set(gca, 'yscale', 'log')
xlabel('t'), ylabel('f(x(t))'), ylim([1e-20, 1])
title('Transformed integrand')
This graph suggests that we capture all of the integrand values that are larger than machine epsilon by integrating in t t t from -4 to 4.
Example 9.7.2 f = lambda x: 1 / (1 + x**2)
x = linspace(-4, 4, 500)
subplot(2, 1, 1)
plot(x, f(x)), yscale('log')
xlabel('x'), ylabel('f(x)'), ylim([1e-16, 1])
title('Original integrand')
xi = lambda t: sinh( pi * sinh(t) / 2 )
dxi_dt = lambda t: pi/2 * cosh(t) * cosh( pi * sinh(t) / 2 )
integrand = lambda t: f(xi(t)) * dxi_dt(t)
subplot(2, 1, 2)
plot(x, integrand(x)), yscale('log')
xlabel('t'), ylabel('f(x(t))'), ylim([1e-16, 1])
title('Transformed integrand')
This graph suggests that we capture all of the integrand values that are larger than machine epsilon by integrating in t t t from -4 to 4.
Function 9.7.1 implements double exponential integration by applying the adaptive integrator Function 5.7.1 to (9.7.7) . It truncates the interval to − M ≤ t ≤ M -M\le t \le M − M ≤ t ≤ M by increasing M M M until the integrand is too small to matter relative to the error tolerance.
Integration over ( − ∞ , ∞ ) (-\infty,\infty) ( − ∞ , ∞ ) About the code
The test isinf(x(M))
in line 17 checks whether x ( M ) x(M) x ( M ) is larger than the maximum double-precision value, causing it to overflow to Inf
.
Integration over ( − ∞ , ∞ ) (-\infty,\infty) ( − ∞ , ∞ ) 1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
function [I, x] = intinf(f, tol)
% INTINF Adaptive doubly exponential integration over (-inf,inf).
% Input:
% f integrand (function)
% tol error tolerance (positive scalar)
% Output:
% I approximation to intergal(f) over (-inf,inf)
% x evaluation nodes (vector)
xi = @(t) sinh(sinh(t));
dxi_dt = @(t) cosh(t) .* cosh(sinh(t));
g = @(t) f(xi(t)) .* dxi_dt(t);
% Find where to truncate the integration interval.
M = 3;
while (abs(g(-M)) > tol/100) || (abs(g(M)) > tol/100)
M = M + 0.5;
if isinf(x(M))
warning("Function may not decay fast enough.")
M = M - 0.5;
break
end
end
[I, t] = intadapt(g, -M, M, tol);
x = xi(t);
end
Integration over ( − ∞ , ∞ ) (-\infty,\infty) ( − ∞ , ∞ ) 1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
def intinf(f, tol):
"""
intinf(f, tol)
Perform doubly-exponential integration of function f over (-Inf,Inf), using
error tolerance tol. Return integral and a vector of the nodes used.
"""
xi = lambda t: np.sinh(np.sinh(t))
dxi_dt = lambda t: np.cosh(t) * np.cosh(np.sinh(t))
g = lambda t: f(xi(t)) * dxi_dt(t)
M = 3
while (abs(g(-M)) > tol/100) or (abs(g(M)) > tol/100):
M += 0.5
if np.isinf(xi(M)):
warnings.warn("Function may not decay fast enough.")
M -= 0.5
break
I, t = intadapt(g,-M,M,tol)
x = xi(t)
return I, x
About the code
The test isinf(x(M))
in line 17 checks whether x ( M ) x(M) x ( M ) is larger than the maximum double-precision value, causing it to overflow to Inf
.
9.7.3 Integrand singularity ¶ If f f f asymptotically approaches infinity as x x x approaches an integration endpoint, its exact integral may or may not be finite. If f f f is integrable, then the part of the integration interval near the singularity needs to be more finely resolved than the rest of it.
Let’s consider
∫ 0 1 f ( x ) d x , \int_0^1 f(x)\,dx, ∫ 0 1 f ( x ) d x , where f f f and/or a derivative of f f f is unbounded at the left endpoint, zero. The change of variable
x ( t ) = 2 1 + exp ( 2 sinh t ) x(t) = \frac{2}{1+\exp(2 \sinh t)} x ( t ) = 1 + exp ( 2 sinh t ) 2 satisfies x ( 0 ) = 1 x(0)=1 x ( 0 ) = 1 and x → 0 + x\to 0^+ x → 0 + as t → ∞ t\to \infty t → ∞ , thereby transforming the integration interval to t ∈ ( 0 , ∞ ) t\in(0,\infty) t ∈ ( 0 , ∞ ) and placing the singularity at infinity. The chain rule implies
∫ 0 1 f ( x ) d x = ∫ 0 ∞ f ( x ( t ) ) d x d t d t = ∫ 0 ∞ f ( x ( t ) ) cosh t cosh ( sinh t ) 2 d t . \begin{split}
\int_{0}^1 f(x)\, dx &= \int_{0}^\infty f(x(t)) \frac{dx}{dt}\, dt \\
&= \int_{0}^\infty f(x(t)) \frac{\cosh t}{\cosh(\sinh t)^2} \, dt.
\end{split} ∫ 0 1 f ( x ) d x = ∫ 0 ∞ f ( x ( t )) d t d x d t = ∫ 0 ∞ f ( x ( t )) cosh ( sinh t ) 2 cosh t d t . Now the growth of f f f and cosh t \cosh t cosh t together are counteracted by the double exponential denominator, allowing easy truncation of (9.7.13) . This variable transformation is paired with adaptive integration in Function 9.7.2 .
Integration with endpoint singularitiesAbout the code
The test iszero(x(M))
in line 17 checks whether x ( M ) x(M) x ( M ) is less than the smallest positive double-precision value, causing it to underflow to zero.
Integration with endpoint singularities1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
function [I, x] = intsing(f, tol)
% INTSING Adaptively integrate a function with a singularity at the left endpoint.
% Input:
% f integrand (function)
% tol error tolerance (positive scalar)
% Output:
% I approximation to integral(f) over (0,1)
% x evaluation nodes (vector)
xi = @(t) 2 ./ (1 + exp( 2*sinh(t) ));
dxi_dt = @(t) cosh(t) ./ cosh( sinh(t) ).^2;
g = @(t) f(xi(t)) .* dxi_dt(t);
% Find where to truncate the integration interval.
M = 3;
while abs(g(M)) > tol/100
M = M + 0.5;
if iszero(x(M))
warning("Function may grow too rapidly.")
M = M - 0.5;
break
end
end
[I, t] = intadapt(g, 0, M, tol);
x = xi(t);
end
Integration with endpoint singularities1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
def intsing(f, tol):
"""
intsing(f, tol)
Adaptively integrate function f over (0,1), where f may be
singular at zero, with error tolerance tol. Returns the
integral estimate and a vector of the nodes used.
"""
xi = lambda t: 2 / (1 + np.exp( 2*np.sinh(t) ))
dxi_dt = lambda t: np.cosh(t) / np.cosh(np.sinh(t))**2
g = lambda t: f(xi(t)) * dxi_dt(t)
# Find where to truncate the integration interval.
M = 3
while abs(g(M)) > tol/100:
M += 0.5
if xi(M) == 0:
warnings.warn("Function may grow too rapidly.")
M -= 0.5
break
I, t = intadapt(g, 0, M, tol)
x = xi(t)
return I, x
About the code
The test iszero(x(M))
in line 17 checks whether x ( M ) x(M) x ( M ) is less than the smallest positive double-precision value, causing it to underflow to zero.
Example 9.7.4 (Singularity at an endpoint)
Let’s use Function 9.7.2 to compute
∫ 0 0.01 1 x d x . \int_0^{0.01} \frac{1}{\sqrt{x}}\, dx. ∫ 0 0.01 x 1 d x . Since the integration interval is not [ 0 , 1 ] [0,1] [ 0 , 1 ] , we must first use the change of variable s = 100 t s=100t s = 100 t , yielding
∫ 0 1 1 10 s d s = 0.2. \int_0^{1} \frac{1}{10\sqrt{s}}\, ds = 0.2. ∫ 0 1 10 s 1 d s = 0.2. In order to use Function 5.7.1 , we must truncate on the left to avoid evaluation at zero, where f f f is infinite. Since the integral from 0 to δ is 20 δ 20\sqrt{\delta} 20 δ , we use δ = ( ϵ / 20 ) 2 \delta=(\epsilon/20)^2 δ = ( ϵ /20 ) 2 to achieve error tolerance ε.
Example 9.7.4 f(x) = 1 / (10 * sqrt(x))
tol = [1 / 10^d for d in 5:0.5:14]
err = zeros(length(tol), 2)
len = zeros(Int, length(tol), 2)
for (i, tol) in enumerate(tol)
I1, x1 = FNC.intadapt(f, (tol/20)^2, 1, tol)
I2, x2 = FNC.intsing(f, tol)
@. err[i, :] = abs(0.2 - [I1, I2])
@. len[i, :] = length([x1, x2])
end
plot(len, err, m=:o, label=["direct" "double exponential"])
n = [30, 3000]
plot!(n, 30n.^(-4);
color=:black, l=:dash,
label="fourth-order", legend=:bottomleft,
xaxis=(:log10, "number of nodes"),
yaxis=(:log10, "error"),
title="Comparison of integration methods")
As in Demo 9.7.3 , the double exponential method is more accurate than direct integration by a few orders of magnitude. Equivalently, the same accuracy can be reached with many fewer nodes.
Example 9.7.4 f = @(x) 1 ./ (10 * sqrt(x));
tol = 1 ./ 10.^(5:0.5:14);
err = zeros(length(tol), 2);
len = zeros(length(tol), 2);
for k = 1:length(tol)
[I1, x1] = intadapt(f, (tol/20)^2, 1, tol);
[I2, x2] = intsing(f, tol);
err(k, :) = abs(0.2 - [I1, I2]);
len(k, :) = [length(x1), length(x2)];
end
clf, loglog(len, err, 'o-')
n = [30, 3000];
hold on, loglog(n, 30 * n.^(-4), 'k--') % 4th order error
legend("direct", "double exponential", "4th order", location="southwest")
title(("Comparison of integration methods"));
As in Demo 9.7.3 , the double exponential method is more accurate than direct integration by a few orders of magnitude. Equivalently, the same accuracy can be reached with many fewer nodes.
Example 9.7.4 f = lambda x: 1 / (10 * sqrt(x))
exact = 0.2
tol = array([1 / 10**d for d in arange(5, 14, 0.5)])
err = zeros((tol.size, 2))
length = zeros((tol.size, 2))
for k in range(tol.size):
I1, x1 = FNC.intadapt(f, (tol[k]/20)**2, 1, tol[k])
I2, x2 = FNC.intsing(f, tol[k])
err[k] = abs(exact - array([I1, I2]))
length[k] = [x1.size, x2.size]
loglog(length, err, "-o")
# plot(len,err,m=:o,label=["direct" "double exponential"])
n = array([100, 10000])
loglog(n, 30 / n**4, 'k--')
xlabel("number of nodes"), ylabel("error")
title("Comparison of integration methods")
legend(["direct", "double exponential", "4th-order"], loc="lower left");
As in Demo 9.7.3 , the double exponential method is more accurate than direct integration by a few orders of magnitude. Equivalently, the same accuracy can be reached with many fewer nodes.
Double exponential integration is an effective general-purpose technique for improper integrals that usually outperforms interval truncation in the original variable. There are specialized methods tailored to specific singularity types that can best it, but those require more analytical work to use properly.
9.7.4 Exercises ¶ ⌨ Use Function 9.7.1 to estimate the given integral with error tolerances 1 0 − 3 , 1 0 − 6 , 1 0 − 9 , 1 0 − 12 10^{-3},10^{-6},10^{-9},10^{-12} 1 0 − 3 , 1 0 − 6 , 1 0 − 9 , 1 0 − 12 . For each result, show the actual error and the number of nodes used.
(a) ∫ − ∞ ∞ 1 1 + x 2 + x 4 d x = π 3 \displaystyle\int_{-\infty}^\infty \dfrac{1}{1+x^2+x^4}\, dx = \dfrac{\pi}{\sqrt{3}} ∫ − ∞ ∞ 1 + x 2 + x 4 1 d x = 3 π
(b) ∫ − ∞ ∞ e − x 2 cos ( x ) d x = e − 1 / 4 π \displaystyle\int_{-\infty}^\infty e^{-x^2}\cos(x)\, dx = e^{-1/4}\sqrt{\pi} ∫ − ∞ ∞ e − x 2 cos ( x ) d x = e − 1/4 π
(c) ∫ − ∞ ∞ ( 1 + x 2 ) − 2 / 3 d x = π Γ ( 1 / 6 ) Γ ( 2 / 3 ) \displaystyle\int_{-\infty}^\infty (1+x^2)^{-2/3}\, dx = \dfrac{\sqrt{\pi}\,\Gamma(1/6)}{\Gamma(2/3)} ∫ − ∞ ∞ ( 1 + x 2 ) − 2/3 d x = Γ ( 2/3 ) π Γ ( 1/6 ) (use gamma()
for Γ ( ) \Gamma() Γ ( ) )
⌨ Use Function 9.7.2 to estimate the given integral, possibly after rewriting the integral into the form (9.7.11) with a left-endpoint singularity. Use error tolerances 1 0 − 3 , 1 0 − 6 , 1 0 − 9 , 1 0 − 12 10^{-3},10^{-6},10^{-9},10^{-12} 1 0 − 3 , 1 0 − 6 , 1 0 − 9 , 1 0 − 12 , and for each result, show the actual error and the number of nodes used.
(a) ∫ 0 1 ( log x ) 2 d x = 2 \displaystyle\int_{0}^1 (\log x)^2\, dx = 2 ∫ 0 1 ( log x ) 2 d x = 2
(b) ∫ 0 π / 4 tan ( x ) d x = π 2 \displaystyle\int_{0}^{\pi/4} \sqrt{\tan(x)}\, dx = \dfrac{\pi}{\sqrt{2}} ∫ 0 π /4 tan ( x ) d x = 2 π
(c) ∫ 0 1 1 1 − x 2 d x = π 2 \displaystyle\int_{0}^1 \frac{1}{\sqrt{1-x^2}}\, dx = \dfrac{\pi}{2} ∫ 0 1 1 − x 2 1 d x = 2 π
For integration on a semi-infinite interval such as x ∈ [ 0 , ∞ ) x\in [0,\infty) x ∈ [ 0 , ∞ ) , another double exponential transformation is useful: x ( t ) = exp ( sinh t ) x(t)=\exp\left( \sinh t \right) x ( t ) = exp ( sinh t ) .
(a) ✍ Show that t ∈ ( − ∞ , ∞ ) t\in(-\infty,\infty) t ∈ ( − ∞ , ∞ ) is mapped to x ∈ ( 0 , ∞ ) x\in (0,\infty) x ∈ ( 0 , ∞ ) .
(b) ✍ Derive an analog of (9.7.7) for the chain rule on ∫ 0 ∞ f ( x ) d x \int_0^\infty f(x)\,dx ∫ 0 ∞ f ( x ) d x .
(c) ✍ Show that truncation of t t t to [ − M , M ] [-M,M] [ − M , M ] will truncate x x x to [ 1 / μ , μ ] [1/\mu,\mu] [ 1/ μ , μ ] for some positive μ.
(d) ⌨ Write a function intsemi(f,tol)
for the semi-infinite integration problem. Test it on the integral
∫ 0 ∞ e − x x d x = π . \displaystyle\int_0^\infty \frac{e^{-x}}{\sqrt{x}}\,dx = \sqrt{\pi}. ∫ 0 ∞ x e − x d x = π .