📂Numerical Analysis

Tricks for Variable Substitution to Compute Improper Integrals Numerically

Theorem ¹

Let’s say $0 < a < b < \infty$.

• [1]: If $ 0 < p < 1$ then $$\int_{0}^{b} {{ f(x) } \over {x^{p} }} dx = \int_{0}^{{{ 1 } \over { 1-p }} b^{1-p} } f \left( \left[ ( 1- p ) m \right]^{{{ 1 } \over { 1-p }}} \right) dm$$
• [2]: If $ 1 < p$ then $$\int_{a}^{ \infty } {{ f(x) } \over {x^{p} }} dx = \int_{0}^{{{ 1 } \over { p-1 }} a^{1-p}}f \left( \left[ ( p-1 ) m \right]^{{{ 1 } \over { 1-p }}} \right) dm$$

Explanation

Improper integrals are always a headache, but you cannot just give up on solving them. The reason integrals are usually hard is that no matter which trick you try to use, those tricks often end up creating new problems. Fortunately, in numerical integration, as long as you have the correct function, you can still perform the calculation without detailed processes. A typical method is substitution, and if there’s a term you can subtract from the denominator, such as $x^{p}$, even if the function is singular, it can be easily solved.

Example

Contour Integration

For example, consider the following improper integral. $$ \int_{0}^{1} {{e^{-x}} \over {\sqrt{x}}} dx \approx 1.493648265625 $$ If we set $\displaystyle m(x) := \int_{0}^{x} {{1} \over {\sqrt{t} }} dt = 2 \sqrt{x}$, then $\displaystyle \sqrt{x} dm = dx$, then $$ \begin{align*} \int_{0}^{1} {{e^{-x}} \over {\sqrt{x}}} dx =& \int_{m=0}^{2} {{e^{-x}} \over {\sqrt{x}}} \sqrt{x} dm \\ =& \int_{0}^{2} e^{-x} dm \\ =& \int_{0}^{2} e^{- m^2 / 4} dm \end{align*} $$ By using such a trick, if we create an integral over a closed interval, we can apply integration techniques. In practice, by using the Trapezoidal Rule, we can calculate the numerical value as follows.

Speed Improvement

On the other hand, as a more interesting example, consider the following case. $$ \int_{0}^{1} e^{-x} \sqrt{x} dx \approx 0.378944691641 $$ The first thought might be ’there’s no need for substitution’. $e^{-x} \sqrt{x}$ does not have singularities in $[0,1]$, so it can be applied as is with the Trapezoidal Rule. However, it’s not differentiable in $x=0$, due to which there’s no error bound, and the convergence speed slows down. Similarly, if we set $\displaystyle m(x) := \int_{0}^{x} {{1} \over {\sqrt{t} }} dt = 2 \sqrt{x}$, then $\displaystyle \sqrt{x} dm = dx$, $$ \begin{align*} \int_{0}^{1} e^{-x} \sqrt{x} dx =& \int_{m=0}^{2} e^{-x} \sqrt{x} \sqrt{x} dm \\ =& \int_{0}^{2} e^{-x} x dm \\ =& \int_{0}^{2} e^{- m^2 / 4} {{m^2} \over {4}} dm \end{align*} $$ As you can see, after performing the substitution integration, there are no singularities in $m=0$. Likewise, by using the Trapezoidal Rule, we can calculate the numeric value as follows. The left is the value calculated without substitution, and the right is the value calculated after substitution.

Even without substitution, convergence can be achieved. However, compared to the terrifying accuracy towards the true value through substitution, the speed of convergence without substitution is noticeably slower. Given the characteristic of valuing speed in numerical analysis, it is beneficial to minimize such losses. When the time comes to choose, one should prefer the meticulous and precise tricks of mathematicians over the simplicity and beauty of the code.

Proof

Strategy: Since the actual calculation will be entrusted to a computer, we just need to convert it into an integral over a closed interval. The methods are essentially the same.

[1]

If we set $$ m := \int_{0}^{x} {{ 1 } \over { t^{p} }} dt = {{ 1 } \over { 1-p }} x^{1-p} $$, then $\displaystyle {{ dm } \over { dx }} = x^{-p} \iff x^{p} dm = dx$, hence $$ \begin{align*} \int_{0}^{b} {{ f(x) } \over {x^{p} }} dx =& \int_{m=0}^{{{ 1 } \over { 1-p }} b^{1-p}} {{ f(x) } \over {x^{p} }} x^{p} dm \\ =& \int_{0}^{{{ 1 } \over { 1-p }} b^{1-p}}f \left( \left[ ( 1- p ) m \right]^{{{ 1 } \over { 1-p }}} \right) dm \end{align*} $$

■

[2]

If we set $$ m := \int_{x}^{\infty} {{ 1 } \over { t^{p} }} dt = - {{ 1 } \over { 1-p }} x^{1-p} $$, then $\displaystyle {{ dm } \over { dx }} = - x^{-p} \iff - x^{p} dm = dx$, hence $$ \begin{align*} \int_{a}^{ \infty } {{ f(x) } \over {x^{p} }} dx =& \int_{m= - {{ 1 } \over { 1-p }} a^{1-p}}^{0} {{ f(x) } \over {x^{p} }} ( - x^{p} ) dm \\ =& \int_{0}^{{{ 1 } \over { p-1 }} a^{1-p}}f \left( \left[ ( p-1 ) m \right]^{{{ 1 } \over { 1-p }}} \right) dm \end{align*} $$

■

Implementation

The following is the Python code that numerically solves the examples.

import numpy as np
 
def f1(x) :
    return np.exp(-x)
 
def f2(x) :
    return np.exp(-x)*np.sqrt(x)
 
def f3(x) :
    return np.exp(-x)*x
 
def I(f,a,b,n) :
    h=(b-a)/n
    x=np.linspace(a,b,n+1)
    x=f(x)
    x[0]=x[0]/2
    x[-1]=x[-1]/2
    return(h*sum(x))
 
def CVI(f,a,b,n) :
    b=2*np.sqrt(b)
    h=(b-a)/n
    x=np.linspace(a,b,n+1)
    x=f((x**2)/4)
    x[0]=x[0]/2
    x[-1]=x[-1]/2
    return(h*sum(x))
 
 
TV1=1.49364826562485405079
TV2=0.37894469164098470380
 
print("")
 
print("True Value 1 : %2.12f" % TV1) 
print("|   n |  In       |  I-In      |")
print("="*32)
for n in range(3,10) :
    print("| %3d |  %8.6f |  %9.2e |" % (
        2**n,
        CVI(f1,0,1,2**n),
        TV1-CVI(f1,0,1,2**n)))
print("")
 
print("True Value 2 : %2.12f" % TV2)
print("|   n |  In       |  I-In      |")
print("="*32)
for n in range(3,10) :
    print("| %3d |  %8.6f |  %9.2e |" % (
        2**n,
        I(f2,0,1,2**n),
        TV2-I(f2,0,1,2**n)))
print("")
 
print("True Value 2 : %2.12f" % TV2)
print("|   n |  In       |  I-In      |")
print("="*32)
for n in range(3,10) :
    print("| %3d |  %8.6f |  %9.2e |" % (
        2**n,
        CVI(f3,0,1,2**n),
        TV2-CVI(f3,0,1,2**n)))
 
print("")