Special techniques - Numerical Analysis

from pylab import *

1Singular integrands¶

Example 1 (Singular integral)

I = \int_0^1 \sqrt{x} \ud x

(1)

The integrand does not have derivatives at $x=0$ .

from scipy.integrate import fixed_quad,trapezoid

f = lambda x: sqrt(x)
a,b = 0.0,1.0
qe = 2.0/3.0 # Exact integral

n,N = 2,10
e1,e2,nodes = zeros(N),zeros(N),zeros(N)
for i in range(N):
    x = linspace(a,b,n)
    val = trapezoid(f(x),dx=(b-a)/(n-1))
    e1[i] = abs(val - qe)
    val = fixed_quad(f,a,b,n=n)
    nodes[i] = n
    e2[i] = abs(val[0]-qe)
    if i>0:
        print('%5d %20.10e %12.4e %20.10e %12.4e' % 
              (n,e1[i],e1[i-1]/e1[i],e2[i],e2[i-1]/e2[i]))
    else:
        print('%5d %20.10e %12.4e %20.10e %12.4e' % 
              (n,e1[i],0,e2[i],0))
    n = 2*n

figure()
loglog(nodes,e1,'o-')
loglog(nodes,e2,'*-')
legend(('Trapezoid','Gauss'))
xlabel('n'), ylabel('Error');

    2     1.6666666667e-01   0.0000e+00     7.2206720124e-03   0.0000e+00
    4     3.5384383294e-02   4.7102e+00     1.1609787082e-03   6.2195e+00
    8     1.0374676944e-02   3.4106e+00     1.6891343351e-04   6.8732e+00
   16     3.3932245922e-03   3.0575e+00     2.2964832819e-05   7.3553e+00
   32     1.1610774091e-03   2.9225e+00     3.0007011894e-06   7.6532e+00
   64     4.0523537954e-04   2.8652e+00     3.8373151734e-07   7.8198e+00
  128     1.4266799811e-04   2.8404e+00     4.8523810858e-08   7.9081e+00
  256     5.0411536176e-05   2.8301e+00     6.1008635921e-09   7.9536e+00
  512     1.7837193144e-05   2.8262e+00     7.6484252265e-10   7.9766e+00

 1024     6.3136745539e-06   2.8252e+00     9.5747521023e-11   7.9881e+00

With trapezoidal rule, we get ( $2^{1.5} \approx 2.825$ )

|E_n(f)| \approx \frac{c}{n^{3/2}}

(2)

while with Gauss quadrature

|E_n(f)| \approx \frac{c}{n^3}

(3)

Gauss quadrature still performs well when the singularity is at the end-point of the interval.

Example 2

\int_0^1 \sqrt{|x - \shalf|} \ud x = \frac{\sqrt{2}}{3} \approx I_n

(6)

f = lambda x: sqrt(abs(x-0.5))
a,b = 0.0,1.0
c = 0.5*(a+b)
qe = sqrt(2.0)/3.0 # Exact integral

n,N = 2,7
err1,err2,nodes = zeros(N),zeros(N),zeros(N)
for i in range(N):
    nodes[i] = n
    val = fixed_quad(f,a,b,n=n)
    err1[i] = abs(val[0]-qe)
    val = fixed_quad(f,a,c,n=n/2) + fixed_quad(f,c,b,n=n/2)
    err2[i] = abs(val[0]+val[2]-qe)
    if i>0:
        print('%5d %14.6e %12.4e %14.6e %12.4e' %
              (n,err1[i],err1[i-1]/err1[i],err2[i],err2[i-1]/err2[i]))
    else:
        print('%5d %14.6e %12.4e %14.6e %12.4e' % (n,err1[i],0,err2[i],0))
    n = 2*n

figure()
loglog(nodes,err1,'o-',nodes,err2,'s-')
xlabel('n'), ylabel('Error')
legend(('Gauss','Gauss+Gauss'));

    2   6.588045e-02   0.0000e+00   2.859548e-02   0.0000e+00
    4   2.572864e-02   2.5606e+00   5.105786e-03   5.6006e+00
    8   9.738642e-03   2.6419e+00   8.209359e-04   6.2195e+00
   16   3.583266e-03   2.7178e+00   1.194398e-04   6.8732e+00
   32   1.294540e-03   2.7680e+00   1.623859e-05   7.3553e+00
   64   4.628656e-04   2.7968e+00   2.121816e-06   7.6532e+00
  128   1.645897e-04   2.8122e+00   2.713392e-07   7.8198e+00

Now Gauss quadrature is not so good because singularity is inside the interval. A trick here is to break it into two intervals and compute each with its own Gauss quadrature

\int_0^\half \sqrt{|x - \shalf|} \ud x + \int_\half^1 \sqrt{x - \shalf} \ud x \approx I_{n/2} + I_{n/2}

(7)

Now the singularity is at the end-points and Gauss quadrature is more accurate.

2Change of variable¶

Many times, a change of variable can improve the regularity of the integrand.

Example 5 (Infinite interval of integration)

I = \int_1^\infty \frac{f(x)}{x^p} \ud x, \qquad p > 1

(12)

Assume $\lim_{x \to \infty} f(x)$ exists and $f(x)$ is smooth on $[1,\infty)$ .

x = \frac{1}{u^\alpha}, \qquad \ud x = -\frac{\alpha}{u^{1+\alpha}} \ud u, \quad \alpha > 0

(13)

I = \alpha \int_0^1 u^{(p-1)\alpha-1} f(1/u^\alpha) \ud u

(14)

Integrand can be singular at $u=0$ , so pick α so that the exponent $(p-1)\alpha-1$ is large $(\ge 1)$ . As a specific example consider

I = \int_1^\infty \frac{f(x)}{x\sqrt{x}} \ud x

(15)

x = \frac{1}{u^4} \quad\implies\quad I = 4 \int_0^1 u f(1/u^4) \ud u

(16)

For $x \to \infty$ , $f(x)$ behaves like

f(x) = c_0 + \frac{c_1}{x} + \frac{c_2}{x^2} + \ldots

(17)

uf(1/u^4) = c_0 u + c_1 u^5 + c_2 u^9 + \ldots

(18)

and the integrand behaves nicely as $u \to 0$ .

Example 6

I = \int_{-1}^1 \frac{f(x)}{\sqrt{1-x^2}} \ud x

(19)

For this weight function, the orthogonal polynomials are Chebyshev polynomials. The integration nodes are roots of $T_n(x) = \cos(n \cos^{-1}x)$

x_j = \cos\left( \frac{2j-1}{2n}\pi \right), \qquad j=1,2,\ldots,n

(20)

and weights

w_j = \frac{\pi}{n}

(21)

so that

I \approx I_n = \frac{\pi}{n} \sum_{j=1}^n f(x_j)

(22)

This corresponds to doing the change of variable $x = \cos\theta$ and applying the mid-point rule in the θ variable

I - \int_0^\pi f(\cos\theta) \ud\theta \approx \frac{\pi}{n} \sum_{j=1}^n f(\cos\theta_j), \qquad \theta_j = \frac{(2j-1)\pi}{2n}

(23)

3Analytic treatment of singularity¶

Consider an integral with a singular term

I = \int_0^b f(x) \log(x) \ud x

(24)

Since $\log(x)$ is singular at $x=0$ , we have to compute this with some care. Divide the interval into two parts

I = \int_0^\epsilon f(x) \log(x) \ud x + \int_{\epsilon}^b f(x) \log(x) \ud x =: I_1 + I_2

(25)

$I_2$ can be computed by some standard technique and $I_1$ is computed semi-analytically. Assume $f(x)$ has convergent Taylor series on $[0,\epsilon]$

f(x) = \sum_{j=0}^\infty a_j x^j, \qquad x \in [0,\epsilon]

(26)

Then

I_1 = \int_0^\epsilon (\sum_j a_j x^j) \log(x) \ud x = \sum_{j=0}^\infty \frac{a_j \epsilon^{j+1}}{j+1} \left[ \log(\epsilon) - \frac{1}{j+1} \right]

(27)

4Product integration¶

Consider an integral with a singular term

I(f) = \int_a^b w(x) f(x) \ud x

(30)

where $w(x)$ is singular and $f(x)$ is a nice function. Let $\{ f_n \}$ be a sequence of approximations to $f$ such that

$|w|$ is integrable.
$\lim_{n\to\infty} \norm{f - f_n}_\infty = 0$
The integrals
$I_n(f) = \int_a^b w(x) f_n(x) \ud x$
(31)
can be analytically computed.

Then

|I(f) - I_n(f)| \le \norm{f-f_n}_\infty \int_a^b |w(x)| \ud x \to 0

(32)

Example 8

I(f) = \int_0^b f(x) \log(x) \ud x

(33)

Let $f_n$ be piecewise linear approximation on a partition

x_j = jh, \quad j=0,1,\ldots,n \qquad\textrm{where}\qquad h = \frac{b}{n}

(34)

and

f_n(x) = \frac{1}{h}[ (x_j - x) f_{j-1} + (x-x_{j-1}) f_j], \qquad x_{j-1} \le x \le x_j

(35)

for which we know the error estimate

\norm{f-f_n}_\infty \le \frac{h^2}{8} \norm{f''}_\infty

(36)

Then

|I(f) - I_n(f)| \le \frac{h^2}{8} \norm{f''}_\infty \int_0^b |\log(x)| \ud x \to 0

(37)

$I_n$ can be computed analytically. In particular, in the first sub-interval $[x_0,x_1]=[0,h]$ we have

\int_{x_0}^{x_1} \log(x) f_n(x) \ud x = \frac{f_0}{h} \int_0^h (h-x)\log(x) \ud x + \frac{f_1}{h} \int_0^h x \log(x) \ud x

(38)

and these are finite since

\int_0^h \log(x) \ud x = h[\log(h)-1], \qquad \int_0^h x\log(x)\ud x = \frac{h^2}{4}[ 2\log(h) - 1 ]

(39)