# Quadrature and Area estimation

Charbel Tannous

► To cite this version:

Charbel Tannous. Quadrature and Area estimation. Master. Analyse Numérique, UBO Brest, France. 2016, pp.23. hal-04980372

HAL Id: hal-04980372

https://hal.science/hal-04980372v1

Submitted on 6 Mar 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Quadrature and Area estimation

C. Tannous

<figure><img src="image_2.png" /></figure>

Université de Brest, Lab-STICC, CNRS-UMR 6285, F-29200 Brest, FRANCE

Area estimation is based on numerical integration (quadrature) through evaluating the integrand at equally spaced points or on a non-uniform grid in order to approach special points of the integrand. Non-uniformly spaced points at which the integrand is evaluated are generally the zeroes of orthogonal polynomials associated with each method (Gauss-Legendre, Gauss-Chebyshev...). Quadrature and associated formulas pertain to four different possibilities: local (simple) or global (composite) rules and open or closed integration interval. Qualitative and quantitative methods are described along with many applications.

PACS numbers: 02.60.Cb, 02.60.-x, 02.40.-k

Keywords: Numerical simulation, solution of equations, Numerical methods (mathematics), Differential geometry

Version: March 6, 2025

## Contents

I. Introduction 2

A. Archimedes and circle area estimation 2

B. Qualitative estimation of integrals 3

C. Solid angle estimation with spherical trigonometry 4

D. Planimeter and Michelson harmonic analyzer 6

II. Quadrature methods 8

A. Trapezoidal Rule 8

B. Composite Simpson Rule 8

C. Gaussian Quadrature 9

D. Gauss-Legendre quadrature 11

E. Gauss-Chebyshev quadrature 11

F. Richardson Extrapolation 11

G. Romberg Integration 13

H. Adaptive Quadrature Methods 13

I. Filon method 14

J. Monte-Carlo quadrature 15

III. Numerical examples 17

A. Estimation of three special integrals 17

B. Additional examples 18

A. Newton-Cotes formulas 19

1. Open and Closed formulas 19

2. Special case $m = 0$ 20

3. Closed Newton-Cotes formulas: $m \neq 0$ 20

4. Open Newton-Cotes formulas: $m \neq 0$ 20

5. Instability at high degree 20

B. Quadrature in the complex plane 21

References 23
# I. INTRODUCTION

Area estimation is very important for assessing property of land owners such as farmers and real estate developers or whether they are geography scientists, mapping engineers or government officials.

Historically [1] a feature of mathematics which has remained fairly constant from the earliest times to the present day is a general view that its aim is to use numbers to solve problems which arise in the world.

Plato defined a particular view of what mathematics was, in that real mathematics is more abstract numbers are no longer numbers of things or measurements of length, area, or time, but have an independent existence as objects which you reason with.

One of Plato longest and clearest mathematical discussions, often referred to, is in the dialogue called the Meno.

<figure><img src="image_3.png" /><figcaption>**Fig.1:** Meno argument related to area estimation is illustrated with a large square with 4 feet side (with $16 \text{ feet}^2$ area), the four smaller squares have each a 2 feet side(4 feet$^2$ area). The four diagonals form a square of $8 \text{ feet}^2$ area.</figcaption></figure>

## A. Archimedes and circle area estimation

Like Plato, Archimedes was interested in numbers and particularly in the evaluation of $\pi$ by estimating it from the area of a circle or perimeter of its circumference.

<figure><img src="image_4.png" /><figcaption>**Fig.2:** Archimedes used the circumference of the triangles with side $r$ and angle $\theta$ inscribed within the circle. The approach can be extended to the circle surface as described by Rice [2].</figcaption></figure>

The question is to find approximations for the circumference $C$ of a circle. One can approximate the value of $C$ with the perimeter of a polygon inscribed within the circle (cf. Fig. 2). Archimedes started with a hexagon inscribed within a circle of radius $r$. The hexagon is made up of six equilateral triangles, each of which has side length $r$. $C \approx r \times 6 = 6r$, and $\pi$ can be obtained from the circumference over the diameter ratio: $\pi \approx 6r/2r = 3$ the Biblical value of $\pi$ in Genesis [1].

In order to improve this approximation for $\pi$, one doubles repetitively the number of sides of the polygon inscribed inside the circle.
Thus an increasing sequence $C_0 = \ell_0 \times 6, C_1 = \ell_1 \times 6 \times 2, C_2 = \ell_2 \times 6 \times 2^2, \dots, C_n = \ell_n \times 6 \times 2^n, \dots$ of approximations for $C$ is obtained. Finally one may define the circumference $C$ by $C = \lim_{n \to \infty} C_n$.

Geometrically, it is possible to use Pythagoras theorem to find the recursion relating the different lengths $\ell_i$ such that:

$$ \ell_{i+1} = r \sqrt{2 - \sqrt{4 - (\ell_i/r)^2}} \quad (1) $$

Iterating 4 times, Archimedes found $3+10/71 < \pi < 3+1/7$, without using square root operations, proceeding rather with geometrical constructions.

Rice [2] described an alternate Archimedes method based on evaluating the surface of triangles in order to estimate $\pi$. The area of each triangle is $\sin(\theta)/2$ and the total area is $2^n \sin(\pi/2^n)$ giving $\lim_{n \to \infty} 2^n \sin(\pi/2^n) = \pi$. During the evaluation, one should not get into a circular argument pitfall of using $\pi$ in a redundant fashion and rather use geometric constructions based on degrees instead of radians in order to estimate $\sin(\pi/2^n)$ via $\sin(180^\circ/2^n)$ as done by Osler [3].

Subdividing into 96 triangles inscribed within the circle, Archimedes [2] finally gets:

$$ 3 \frac{1137}{8069} < \pi < 3 \frac{1335}{9347} \quad \text{or} \quad 3.1409 < \pi < 3.1428 \quad (2) $$

## B. Qualitative estimation of integrals

Migdal [4] offers several methods for estimating integrals. He starts with the area of a right triangle by supposing that it must be of the form $a^2 f(\alpha)$ where $a$ is the hypotenuse and $f$ a scaling function of the apex angle $\alpha$. Then he obtains from decomposing the triangle Pythagoras theorem since the total area is $a^2 f(\alpha) = b^2 f(\alpha) + c^2 f(\alpha)$ yielding $a^2 = b^2 + c^2$

In the arbitrary triangle case, its area (cf. Fig. 3) is given by: $bh/2 = c^2 f(\alpha) + a^2 f(\gamma)$ where $b, h$ are base and height of the triangle.

<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /><figcaption>**Fig.3:** (Left) Right triangle decomposed into smaller homothetic right triangles. (Right) Arbitrary triangle decomposed into smaller right triangles</figcaption></figure>

Approximating an integral is based on simplifying the integrand, expanding it, examining its singularities... Abramowitz [5] has given the following indications [6]:

* Confirm the existence of the integral,

* Ascertain the important ranges of the parameters involved,

* Reduce the integral to its simplest form,

* Determine the essential parameters which are involved,

* Determine the accuracy to which numerical values (if desired) are to be given.

An example is the BCS integral [7] of superconductivity giving the critical temperature $T_c$:

$$ \int_{0}^{\beta_c \hbar \omega_c / 2} \frac{\tanh(x)}{x} dx, \quad \beta_c = 1/k_B T_c \quad (3) $$

where $k_B$ is Boltzmann constant and $\omega_c$ a characteristic pulsation [7]. A very lengthy integration [8] yields $\ln(a\beta_c\hbar\omega_c)$ with $a = 2e^\gamma/\pi \approx 1.13$ where $\gamma = 0.577...$ is Euler-Mascheroni [5] constant given by: $\lim_{p \to \infty} 1 + \frac{1}{2} + \frac{1}{3} + \frac{1}{4} + \dots + \frac{1}{p} - \ln p$.
An estimation of this integral can be made [9] in an easy way by simplifying the tanh(x) function as a linear function over the interval [0, A] with A a cutoff above which it is one i.e.

$$ \tanh(x) = \begin{cases} x, & x \in [0, A[; \\ 1, & x \ge A \end{cases} \qquad (4) $$

then the integral is written as:

$$ \int_{0}^{\beta_c \hbar \omega_c / 2} \frac{\tanh(x)}{x} dx \sim \int_{0}^{A} dx + \int_{A}^{\beta_c \hbar \omega_c / 2} \frac{1}{x} dx = A + \ln(\beta_c \hbar \omega_c / 2A) \qquad (5) $$

yielding the crucial logarithmic behavior.

### C. Solid angle estimation with spherical trigonometry

Another example is the estimation of a solid angle given by the area on a unit sphere as displayed in Fig. 4:

<figure><img src="image_7.png" /><figcaption>Fig. 4: The solid angle depicted on the unit sphere is $d\Omega = dA/r^2$ where the area element $dA = r^2 \sin\theta d\theta d\phi$ is on the radius r sphere.</figcaption></figure>

A solid angle is an area on the unit sphere is given by (cf. Fig. 4):

$$ \Omega = \int_{\phi_{min}}^{\phi_{max}} d\phi \int_{\theta_{min}}^{\theta_{max}} d\theta \sin\theta \qquad (6) $$

where $\phi_{min}, \phi_{max}, \theta_{min}, \theta_{max}$ are the limiting values of the respective $\phi, \theta(\phi)$ angles that are inter-dependent.

Spherical trigonometry [10] is extremely helpful in evaluation spherical areas without going through tedious 2D integration on a sphere. Casey [11] is a fine book of spherical trigonometry that helps to find the solid angle subtended by a spherical triangle (cf. Fig. 5) as:

$$ \cos\left(\frac{1}{2}\Omega\right) = \frac{1 + \cos a + \cos b + \cos c}{4 \cos \frac{1}{2} a \cos \frac{1}{2} b \cos \frac{1}{2} c} \qquad (7) $$

where a, b, c are the arc lengths of the triangle ABC (cf. eq. (359) in Casey [11] book).

It is possible to bypass arc length evaluation and express $\Omega$ simply with vectors [12] $R_1, R_2, R_3$ originating from sphere center O (cf. Fig. 5) and linking each apex of the plane triangle $R_1R_2R_3$ to its corresponding one belonging to the spherical triangle ABC:

$$ \tan\left(\frac{1}{2}\Omega\right) = \frac{[R_1, R_2, R_3]}{R_1 R_2 R_3 + (R_1 \cdot R_2) R_3 + (R_1 \cdot R_3) R_2 + (R_2 \cdot R_3) R_1} \qquad (8) $$

where $[R_1, R_2, R_3]$ is a triple product $R_1 \cdot (R_2 \times R_3)$.

An interesting application of the triangular solid angle is the extension to one subtended by a rectangle, i.e. a screen that belongs to one of four classes:
<figure><img src="image_8.png" /><figcaption>**Fig.5:** Spherical triangle ABC on the unit sphere giving the corresponding solid angle subtending the plane triangle R₁R₂R₃. Spherical triangle ABC is displayed as an inset with its arc lengths a, b, c.</figcaption></figure>

1. Desktop computer Monitor

2. Home Television screen

3. Cinema or Wall (building facades) screen

4. Hand-held, augmented/Virtual-Reality (gears, helmets...) and 3D displays

In this case, one considers a rectangular pyramid solid angle extending formula eq. 8 to a rectangle made of four sets of two right triangles with sides a, b (displayed within a blue rectangle in Fig. 6), thus one obtains:

$$ \Omega = 4 \tan^{-1} \frac{ab}{h\sqrt{a^2 + b^2 + h^2}} \quad (9) $$

This formula satisfies the limits: $lim_{h→0} Ω = 2π$, $lim_{h→∞} Ω = 0$, $lim_{a,b→0} Ω = 0$ and $lim_{a,b→∞} Ω = 2π$. It is fully verified by calculating the solid angle with Monte-Carlo quadrature (cf. Section II J) as displayed in Fig. 7.

<figure><img src="image_9.png" /><figcaption>**Fig.6:** Rectangular screen subtending a solid angle from point O at a viewing distance h perpendicularly to the rectangle at its center. Right triangles having a, b sides are shown in the blue rectangle.</figcaption></figure>

Fig. 7 shows a strong disagreement between Monte-Carlo quadrature and the analytical result $Ω = 4 \tan^{-1} \frac{ab}{h\sqrt{a^2+b^2+h^2}}$ over the full interval $h \in [0, 25]$. It is possible to resolve this issue by several methods [13, 14] pertaining to variance reduction methods such as importance sampling, stratified sampling, function approximation, VEGAS algorithm [14] an improved variant of Monte-Carlo method...

The origin of this issue can be revealed by examining the surface integral eq. 6 obtained from the four sets of two right triangles with sides a, b (cf. Fig. 6):

$$ \Omega = 4 \left[ \int_{0}^{\tan^{-1}(b/a)} d\phi \int_{0}^{\cos^{-1}\frac{h}{\sqrt{h^2+a^2/\cos^2\phi}}} d\theta \sin\theta + \int_{\tan^{-1}(b/a)}^{\pi/2} d\phi \int_{0}^{\cos^{-1}\frac{h}{\sqrt{h^2+b^2/\sin^2\phi}}} d\theta \sin\theta \right] \quad (10) $$

Performing integration with respect to $θ$ yields:

$$ \Omega = 4 \left[ \int_{0}^{\tan^{-1}(b/a)} \left( 1 - \frac{h}{\sqrt{h^2 + a^2 / \cos^2 \phi}} \right) d\phi + \int_{\tan^{-1}(b/a)}^{\pi/2} \left( 1 - \frac{h}{\sqrt{h^2 + b^2 / \sin^2 \phi}} \right) d\phi \right] \quad (11) $$
<figure><img src="image_10.png" /><figcaption>**Fig.7:** Solid angle Ω subtended by a rectangle versus observation distance h evaluated by Monte-Carlo quadrature of the surface integral given by eq. 10 and compared to the exact formula $4\tan^{-1}\frac{ab}{h\sqrt{a^2+b^2+h^2}}$. a = 2, b = 1 and the number of Monte-Carlo steps is 10⁵ yielding an average error that is large over the full interval $h \in [0, 25]$. This discrepancy can be treated with variance reduction methods (cf. Section II J).</figcaption></figure>

that can be transformed into:

$$ \Omega = 4 \left[ \pi/2 - \int_{0}^{\tan^{-1}(b/a)} \frac{h}{\sqrt{h^2 + a^2/\cos^2 \phi}} d\phi - \int_{\tan^{-1}(b/a)}^{\pi/2} \frac{h}{\sqrt{h^2 + b^2/\sin^2 \phi}} d\phi \right] \quad (12) $$

Thus it is possible to consider a better approximation such as eq. 12 to represent the surface integral eq. 10. Performing Monte-Carlo quadrature using this approximation leads to a better approximation of Ω over the full interval $h \in [0, 25]$ as displayed in Fig. 8.

<figure><img src="image_11.png" /><figcaption>**Fig.8:** Perfect agreement, for all distances h, between Monte-Carlo quadrature of functions given in eq. 12 and the analytical result $\Omega = 4\tan^{-1}\frac{ab}{h\sqrt{a^2+b^2+h^2}}$. a = 2, b = 1 and the number of Monte-Carlo steps is 10⁵.</figcaption></figure>

## D. Planimeter and Michelson harmonic analyzer

A Mechanical Integrator is a mechanical device that can be used to obtain the continuous quadrature of a curve representing a function with respect to an input variable.

Planimeters and Harmonic analyzers were developed to evaluate Fourier coefficients from periodic functions [15].
Various methods have been proposed for calculating the coefficients of a Fourier series representation of a function that is given in graphical form.

One method is graphical-mechanical, the ordinates being scanned and the area of the curve so obtained being evaluated with a planimeter. Mader [16] has developed a harmonic analyzer that is purely mechanical and that in its most complicated (and expensive) form gives a total of 14 harmonic components with one tracking of the plotted function curve giving an approximate value of a sine or cosine coefficient directly.

Values of the higher order harmonics are obtained in the same way and can be used as correction terms to improve the accuracy of the result.

<figure><img src="image_12.png" /><figcaption>Fig.9: W.F. Stanley apparatus used in mechanical integration to extract Fourier coefficients directly from a periodic curve.</figcaption></figure>

The operation of a linear planimeter can be understood by applying Green theorem [17].

The connection with Green theorem [18] in the plane is formulated as:

Let $P(x, y)$ and $Q(x, y)$ be continuous and have continuous partial derivatives in a region $\mathcal{R}$ and on its boundary $C$. Green theorem states that for a vector field $V = (P, Q)$, its circulation along a curve $C$ is given by:

$$ \oint_C \mathbf{V} \cdot d\mathbf{r} = \oint_C P dx + Q dy = \iint_{\mathcal{R}} \left( \frac{\partial Q}{\partial x} - \frac{\partial P}{\partial y} \right) dxdy \quad (13) $$

The theorem is valid for both simply and multiply connected regions (cf. Fig. 10).

<figure><img src="image_13.png" /><figcaption>Fig.10: At left, singly connected region $\mathcal{R}$ bounded by curve $C$ and at right doubly connected region $\mathcal{R}$ bounded by two curves $C_1$ and $C_2$.</figcaption></figure>

History of harmonic analyzers dates back to Lord Kelvin, Henrici, Sharp, Yule, Michelson and Stratton [19].

A visual detailed book that celebrates a harmonic analyzer designed in the late nineteenth century by physicist Albert Michelson is available [20]. A harmonic analyzer can carry out two related tasks: it can add together weighted sines or cosines to produce a function, and it can perform the inverse operation of decomposing a given function into its constituent sinusoids. This operation is in fact Fourier synthesis.

Michelson and Stratton [19] describe an 80-element harmonic analyzer, a machine with four times as many elements as the analyzer shown previously in the book [20]. They outline the essential mechanical elements of the analyzer, show pages of sample output, and take a brief look at the mathematical approximations and simplifications underlying the machine operation.
## II. QUADRATURE METHODS

Given a function $f(x)$ defined over $[a, b]$ and known at equally spaced points $x_i$, for $i = 0, \dots, n$, where $x_0 = a$ and $x_n = b$, quadrature is considered as local (simple) if one is interested in a given interval such as $[a, a+nh]$ where $h$ is the integration step (usually small since typically $h = \frac{b-a}{n}$). When the full [21] $[a, b]$ interval is considered for integration, quadrature is qualified as global (composite).

Other qualifiers are open when the end points $a, b$ are not considered and closed when $a, b$ are included in the integration process.

### A. Trapezoidal Rule

Suppose that $f : [a, b] \to \mathbf{R}$ is a twice continuously differentiable [21] function on $[a, b]$. Then

$$ \int_{a}^{b} f(x) dx = \frac{f(a) + f(b)}{2} (b - a) - \frac{f''(\xi) (b - a)^3}{12}, \quad \xi \in [a, b] $$

1. The trapezoidal rule is a closed Newton-Cotes formula (cf. Appendix).

2. If $|a - b|$ is small, $\int_{a}^{b} f(x) dx \approx \frac{f(a)+f(b)}{2} (b - a)$ with a truncation error $-f''(\xi) \frac{(b-a)^3}{2}$.

**Theorem 1 (Simpson Rule).** Suppose that $f : [a, b] \to \mathbf{R}$ is a four times continuously differentiable function on $[a, b]$. Then

$$ \int_{a}^{b} f(x) dx = \frac{f(a) + 4f\left(\frac{a+b}{2}\right) + f(b)}{6} (b - a) - \frac{f^{(4)}(\xi) (b - a)^{5}}{2880}, \quad \xi \in [a, b] $$

1. The Simpson rule [21] is also a closed Newton-Cotes formula (cf. Appendix).

2. If $|a - b|$ is small, $\int_{a}^{b} f(x) dx \approx \frac{1}{2} (f(a) + f(\frac{a+b}{2}) + f(b)) (b - a)$ and the truncation error is $-\frac{f^{(4)}(\xi)}{2880} (b - a)^5$.

**Theorem 2 (Midpoint Rule).** Suppose that $f : [a, b] \to \mathbf{R}$ is twice continuously differentiable function on $[a, b]$. Then

$$ \int_{a}^{b} f(x) dx = f\left(\frac{a+b}{2}\right) (b-a) + \frac{f''(\xi)}{24} (b-a)^3, \quad \xi \in [a, b] $$

1. The midpoint rule is an open Newton-Cotes formula (cf. Appendix).

2. If $|a - b|$ is small, $\int_{a}^{b} f(x) dx \approx f(\frac{a+b}{2}) (b - a)$ and the truncation error is $\frac{f''(\xi)}{24} (b - a)^3$.

### B. Composite Simpson Rule

**Theorem 3 (Composite Simpson Rule).** Let $f : [a, b] \to \mathbf{R}$ be a four times continuously differentiable function [21]. Let $n = 2m$, $h = (b - a)/n$ and $x_j = a + jh$ for $j = 0, 1, \dots, n$. then

$$ \int_{a}^{b} f(x) dx = \frac{h}{3} \left( f(x_0) + 2 \sum_{j=1}^{m-1} f(x_{2j}) + 4 \sum_{j=0}^{m-1} f(x_{2j+1}) + f(x_n) \right) - \frac{f^{(4)}(\xi) (b-a)}{180} h^4, \quad \xi \in [a, b] $$

Given: $\int_{a}^{b} f(x) dx \approx \frac{h}{3} \left( f(x_0) + 2 \sum_{j=1}^{m-1} f(x_{2j}) + 4 \sum_{j=0}^{m-1} f(x_{2j+1}) + f(x_n) \right)$ with the truncation error $-\frac{f^{(4)}(\xi) (b-a)}{180} h^4$.

**Example 1:** Use the composite Simpson rule [21] to approximate $\int_{0}^{1} e^{x^2} dx$. Choose the number of subintervals such that the magnitude of the truncation error is smaller than $10^{-4}$.
First choose $m$ such that the magnitude of the truncation error in the composite Simpson rule (Theorem 3) is smaller than $10^{-4}$; namely,

$$ \left| \frac{f^{(4)}(\xi)(b-a)}{180} h^4 \right| < 10^{-4}. $$

Thus $f(x) = e^{x^2}$, $a = 0$, $b = 1$, $x_j = jh$ and $h = 1/n$ where $n = 2m$. Since $|f^{(4)}(x)| = 4e^{x^2}(3+12x^2+4x^4)$, $|f^{(4)}(x)| \le 76e$ on $[0, 1]$. Thus, the magnitude of the truncation error is at most $\frac{76e}{180(2m)^4} = \frac{19e}{720m^4}$. Choose $m$ such that $\frac{19e}{720m^4} < 10^{-4}$; namely, $m > 5.175220955...$. With $m = 6$, $h = 1/12$ and

$$ \begin{aligned} \int_0^1 e^{x^2} dx &\approx \frac{1}{36} \left( f(0) + 2 \sum_{j=1}^5 f(j/6) + 4 \sum_{j=0}^5 f((2j+1)/12) + f(1) \right) \\ &= \frac{1}{36} \left( 1 + 2 \left( e^{(2/12)^2} + \dots + e^{(10/12)^2} \right) + 4 \left( e^{(1/12)^2} + \dots + e^{(11/12)^2} \right) + e \right) \approx 1.46267. \end{aligned} $$

### C. Gaussian Quadrature

Let $f: ]a, b[ \to \mathbb{R}$ be an integrable function on $]a, b[$. $f$ may not be a nice function to integrate. For instance, $f$ may not be bounded at the endpoints, thus $\int_a^b f(x) dx$ is an improper integral [21].

Assuming it is possible to write $f$ as the product of two functions $g$ and $w$, where $g$ is a nice function on $]a, b[$ and $w$ is a function on $]a, b[$ taking only non-negative values (and almost everywhere non-null).

A Gaussian quadrature is a formula of the form

$$ \int_a^b f(x) dx = \int_a^b g(x) w(x) dx \approx \sum_{i=1}^n c_i g(x_i), \qquad (14) $$

where choose nodes $x_1, x_2, \dots, x_n$ are chosen in $[a, b]$ and the weights $c_1, c_2, \dots, c_n$ such that the formula is exact for all polynomials $g$ of degree less than a given constant $k$ (usually $k = 2n$).

It is easy to find a Gaussian quadrature that is exact for polynomials of degree less than $n$. Given any nodes $a \le x_1 < x_2 < \dots < x_n \le b$, let

$$ \ell_i(x) = \prod_{\substack{j=1 \\ j \neq i}}^n \frac{x-x_j}{x_i-x_j} \quad \text{and} \quad c_i = \int_a^b \ell_i(x) w(x) dx \qquad (15) $$

for $1 \le i \le n$. Since any polynomial $p$ of degree less than $n$ can be written as

$$ p(x) = \sum_{i=1}^{n} p(x_i) \ell_i(x), $$

thus:

$$ \int_a^b p(x)w(x) dx = \sum_{i=1}^n p(x_i) \left( \int_a^b \ell_i(x)w(x) dx \right) = \sum_{i=1}^n c_i p(x_i). $$

### Gaussian nodes and weights

Numerical Recipes [14] and Abramowitz-Stegun [5] give the following:

**Example 2:** Find $x_1$ and $c_1$ such that

$$ \int_a^b f(x) dx \approx c_1 f(x_1) $$
<table><caption>**Table 1:** Gaussian nodes and weights of Gaussian integration. $P_n(x)$ are Legendre polynomials [5] whose zeroes are the nodes.</caption><thead><tr><th>n</th><th>Nodes x<sub>j</sub></th><th>weights c<sub>j</sub></th></tr></thead><tbody><tr><td>2</td><td>± 0.5773502692</td><td>1.0</td></tr><tr><td>3</td><td>± 0.7745966692</td><td>0.5555555556</td></tr><tr><td></td><td></td><td>0.0 0.8888888889</td></tr><tr><td>4</td><td>± 0.8611363116</td><td>0.3478548451</td></tr><tr><td></td><td>± 0.3399810436</td><td>0.6521451549</td></tr><tr><td>5</td><td>± 0.9061798459</td><td>0.2369268851</td></tr><tr><td></td><td>± 0.5384693101</td><td>0.4786286705</td></tr><tr><td></td><td></td><td>0.0 0.5688888889</td></tr></tbody></table>

is exact for polynomial of degree less than 2. The node $x_1$ and the weight $c_1$ must satisfy $\int_a^b 1 dx = c_1$ and $\int_a^b x dx = c_1x_1$; namely, $a+b=c_1$ and $\frac{b^2}{2} - \frac{a^2}{2} = c_1x_1$. The values of $c_1$ and $x_1$ satisfying these two equations are $c_1 = b-a$ and $x_1 = (b+a)/2$. The quadrature formula is therefore

$$\int_a^b f(x) dx \approx (b-a)f(\frac{b+a}{2}).$$

This is the Midpoint rule. The truncation error in $f(x) = x^2$ case is $-f''(\xi)(b-a)^3/24$ since:

$$\int_a^b x^2 dx - (b-a) \left(\frac{b+a}{2}\right)^2 = \frac{b^3 - a^3}{3} - \frac{(b-a)(b+a)^2}{4} = -\frac{(b-a)^3}{12}.$$

**Example 3:** Find $x_1, x_2, c_1$ and $c_2$ such that

$$\int_a^b f(x) dx \approx c_1 f(x_1) + c_2 f(x_2)$$

is exact for polynomial of degree less than 4. The nodes $x_1, x_2$ and the weights $c_1, c_2$ must satisfy $\int_a^b 1 dx = c_1 + c_2$, $\int_a^b x dx = c_1x_1+c_2x_2$, $\int_a^b x^2 dx = c_1x_1^2+c_2x_2^2$ and $\int_a^b x^3 dx = c_1x_1^3+c_2x_2^3$; Namely, $b-a = c_1+c_2$, $\frac{b^2-a^2}{2} = c_1x_1+c_2x_2$, $\frac{b^3-a^3}{3} = c_1x_1^2+c_2x_2^2$ and $\frac{b^4-a^4}{4} = c_1x_1^3+c_2x_2^3$. The values of $c_1, c_2, x_1$ and $x_2$ satisfying these four nonlinear equations are $c_1 = c_2 = \frac{b-a}{2}$, $x_1 = z$ and $x_2 = a+b-z$, where $z$ is the positive root of $6z^2 - 6(a+b)z + a^2 + b^2 + 4ab$.

If $a = -1$ and $b = 1$, $c_1 = c_2 = 1$ and $x_1 = -x_2 = 1/\sqrt{3}$. Hence,

$$\int_{-1}^{1} f(x) dx \approx f\left(\frac{1}{\sqrt{3}}\right) + f\left(-\frac{1}{\sqrt{3}}\right).$$

In fact this is the Gauss-Legendre quadrature formula for $n=2$.

In general, choosing nodes $x_1, x_2, \dots, x_n$ and weights $c_1, c_2, \dots, c_n$ such that (14) is exact for polynomials of degree less than $2n$. This is proven in the following.

**Theorem 4.** Let $\{P_0, P_1, P_2, \dots\}$ be an orthogonal set of polynomials on $[a, b]$ with respect to a weight function $w$. Suppose that $P_n$ is of degree exactly $n$. If $p$ is a polynomial of degree less than $2n$, then

$$\int_a^b p(x)w(x) dx = \sum_{j=1}^n c_j p(x_j), \quad c_j = \int_a^b \left( \prod_{\substack{i=1 \\ i \neq j}}^n \frac{x-x_i}{x_j-x_i} \right) w(x) dx$$

and $x_1, x_2, \dots, x_n$ are the roots of the Legendre polynomial $P_n$ of degree $n$.

Substituting $g(x) = 1$ in a Gaussian quadrature formula

$$\int_a^b g(x)w(x) dx \approx \sum_{j=1}^n c_j g(x_j)$$

which is exact for polynomials of degree less than $2n$, leads to: $\int_a^b w(x) dx = \sum_{j=1}^n c_j$.
### D. Gauss-Legendre quadrature

Considering the Legendre polynomials and $w(x) = 1$ for $-1 \le x \le 1$ in Theorem 4, yields the Gauss-Legendre quadrature.

The integral $\int_{-1}^{1} f(x) dx$ is approximately equal to $\sum_{j=1}^{n} c_j f(x_j)$, where the $x_j$'s are the roots of the Legendre polynomial $P_n$ and the $c_j$'s are given by Theorem 4 with $a = -1$, $b = 1$ and $w(x) = 1$ for $-1 \le x \le 1$. The values of $c_j$ and $x_j$ for $n = 2, 3, 4$ and $5$ are given in the appendix.

**Example 4:** Use Gauss-Legendre quadrature with $n = 3$ to approximate

$$ \int_{1}^{3} \frac{\sin^2(x)}{x} dx . $$

Using the change of variable $t = x - 2$, gives:

$$ \begin{aligned} \int_{1}^{3} \frac{\sin^2(x)}{x} dx &= \int_{-1}^{1} \frac{\sin^2(t+2)}{t+2} dt \approx \sum_{i=1}^{3} c_i \frac{\sin^2(x_i+2)}{x_i+2} \\ &\approx 0.5555555556 \frac{\sin^2(2+0.7745966692)}{2+0.7745966692} + 0.8888888889 \frac{\sin^2(2)}{2} \\ &\quad + 0.5555555556 \frac{\sin^2(2-0.7745966692)}{2-0.7745966692} \approx 0.79465267. \end{aligned} $$

In general, to transform an integral of the form $\int_{a}^{b} f(x) dx$ into an integral of the form $\int_{-1}^{1} g(t) dt$, one uses the substitution $t = \frac{x-(a+b)/2}{(b-a)/2}$, where $(a+b)/2$ is the middle of the interval $[a, b]$ and $(b-a)/2$ is half the length of $[a, b]$.

### E. Gauss-Chebyshev quadrature

Using the Chebyshev polynomials and $w(x) = 1/\sqrt{1-x^2}$ for $-1 < x < 1$ in Theorem 4, leads to Gauss-Chebyshev quadrature. The integral $\int_{-1}^{1} \frac{g(x)}{\sqrt{1-x^2}} dx$ is approximately equal to $\sum_{j=1}^{n} c_j g(x_j)$, where the $x_j$'s are the roots of the Chebyshev polynomial $T_n$ and the $c_j$'s are given by Theorem 4 with $a = -1$, $b = 1$ and $w(x) = 1/\sqrt{1-x^2}$ for $-1 < x < 1$. So, $x_j = \cos((2j-1)\pi/(2n))$ for $1 \le j \le n$, and one can prove that $c_j = \pi/n$ for all $j$.

**Example 5:** Use Gauss-Chebyshev quadrature to approximate

$$ \int_{-1}^{1} \frac{x^2}{\sqrt{1-x^2}} dx . $$

Considering Gauss-Chebyshev quadrature with $n = 2$ gives:

$$ \int_{-1}^{1} \frac{x^2}{\sqrt{1-x^2}} dx = c_1 x_1^2 + c_2 x_2^2 = \frac{\pi}{2} (\cos(\pi/4))^2 + \frac{\pi}{2} (\cos(3\pi/4))^2 = 1.57079632679 \dots $$

because Gauss-Chebyshev quadrature with $n = 2$ is exact for polynomial of degree less than $2n = 4$. In general, $\int_{-1}^{1} \frac{p(x)}{\sqrt{1-x^2}} dx = c_1 p(x_1) + c_2 p(x_2)$ for any polynomial $p(x)$ of degree less than 4.

### F. Richardson Extrapolation

Richardson extrapolation is also called extrapolation to the limit [21].

Let $f: \mathbf{R} \to \mathbf{R}$ be a sufficiently continuously differentiable function and $L(f) = f'(c)$ for some $c \in \mathbf{R}$ (or $L(f) = \int_{a}^{b} f(x) dx$). Suppose that $L_h(f)$ is an approximation of $L(f)$ satisfying

$$ L(f) = L_h(f) + \sum_{j=1}^{\infty} K_j h^{2j} \qquad (16) $$
for h near the origin. A procedure that generates better approximations of L(f) than L<sub>h</sub>(f) from a truncation point of view is described below [21].

Replacing h in (16) by h/2 and h/2<sup>2</sup>, yields respectively:

$$
L(f) = L_{h/2}(f) + \sum_{j=1}^{\infty} K_j \left(\frac{h}{2}\right)^{2j} . \qquad (17)
$$

and

$$
L(f) = L_{h/2^2}(f) + \sum_{j=1}^{\infty} K_j \left(\frac{h}{2^2}\right)^{2j} . \qquad (18)
$$

Subtracting (16) from 4 × (17) and dividing the result by 4-1 i.e. 3, leads to:

$$
\begin{align*}
L(f) &= L_{h/2}^{1}(f) + \frac{1}{3} \left( 4 \sum_{j=1}^{\infty} K_j \left(\frac{h}{2}\right)^{2j} - \sum_{j=1}^{\infty} K_j h^{2j} \right) \\
&= L_{h/2}^{1}(f) + \frac{1}{3} \left( \sum_{j=1}^{\infty} \left( 4K_j \left(\frac{h}{2}\right)^{2j} - 4^j K_j \left(\frac{h}{2}\right)^{2j} \right) \right) \\
&= L_{h/2}^{1}(f) + \sum_{j=2}^{\infty} \frac{4-4^j}{3} K_j \left(\frac{h}{2}\right)^{2j}, \tag{19}
\end{align*}
$$

where

$$
L_{h/2}^{1}(f) = \frac{4L_{h/2}(f) - L_{h}(f)}{3}
$$

is an approximation of L(f) with truncation error −K₂h⁴/4 + O(h⁶). In fact, the expression O(hk) replaces a function g(h) for which there exists a constant M such that |g(h)| ≤ Mhk for h close to the origin. L₁ʰ/₂(f) is a better approximation of L(f) than Lʰ/₂(f) since, for h < 1, the truncation error of L₁ʰ/₂(f) in (19) is generally smaller than the one pertaining to Lʰ/₂(f) in (17).

Subtracting (17) from 4 × (18) and dividing the result by 3, yields:

$$
\begin{align*}
L(f) &= L_{h/2^2}^{1}(f) + \frac{1}{3} \left( 4 \sum_{j=1}^{\infty} K_j \left(\frac{h}{2^2}\right)^{2j} - \sum_{j=1}^{\infty} K_j \left(\frac{h}{2}\right)^{2j} \right) \\
&= L_{h/2^2}^{1}(f) + \frac{1}{3} \left( \sum_{j=1}^{\infty} \left( 4K_j \left(\frac{h}{2^2}\right)^{2j} - 4^j K_j \left(\frac{h}{2^2}\right)^{2j} \right) \right) \\
&= L_{h/2^2}^{1}(f) + \sum_{j=2}^{\infty} \frac{4-4^j}{3} K_j \left(\frac{h}{2^2}\right)^{2j}, \tag{20}
\end{align*}
$$

where

$$
L_{h/2^2}^{1}(f) = \frac{4L_{h/2^2}(f) - L_{h/2}(f)}{3}
$$

is an approximation of L(f) with truncation error −K₂h⁴/4³ + O(h⁶).

It can be proven by induction [21] that

$$
L(f) = L_{h/2^k}^1(f) + \sum_{j=2}^{\infty} \frac{4-4^j}{3} K_j \left(\frac{h}{2^k}\right)^{2j}, \quad (21)
$$

where

$$
L_{h/2^k}^{1}(f) = \frac{4L_{h/2^k}(f) - L_{h/2^{k-1}}(f)}{3}
$$
is an approximation of $L(f)$ with truncation error $O(h^4)$.

Subtracting (19) from $4^2 \times (20)$ and dividing the result by $4^2 - 1$ i.e. 15, yields:

$$ 
\begin{aligned}
L(f) &= L_{h/2^2}^2(f) + \frac{1}{15} \left( 4^2 \sum_{j=2}^{\infty} \frac{4-4^j}{3} K_j \left( \frac{h}{2^2} \right)^{2j} - \sum_{j=2}^{\infty} \frac{4-4^j}{3} K_j \left( \frac{h}{2} \right)^{2j} \right) \\
&= L_{h/2^2}^2(f) + \frac{1}{15} \sum_{j=2}^{\infty} \left( 4^2 \frac{4-4^j}{3} K_j \left( \frac{h}{2^2} \right)^{2j} - 4^j \frac{4-4^j}{3} K_j \left( \frac{h}{2^2} \right)^{2j} \right) \\
&= L_{h/2^2}^2(f) + \sum_{j=3}^{\infty} \frac{(4^2-4^j)(4-4^j)}{(15)(3)} K_j \left( \frac{h}{2^2} \right)^{2j},
\end{aligned} 
\quad (22) $$

where

$$ L_{h/2^2}^2(f) = \frac{4^2 L_{h/2^2}^1(f) - L_{h/2}^1(f)}{15} $$

is an approximation [21] of $L(f)$ with truncation error $K_3h^6/64 + O(h^8)$, assuming that $L_{h/2^2}^2(f)$ is the best approximation of $L_h(f)$ obtained. For small $h < 1$, the truncation error is generally smaller for $L_{h/2^2}^2(f)$ than other approaches.

### G. Romberg Integration

Romberg integration is nothing else than Richardson extrapolation with $L(f) = \int_a^b f(x) dx$ and

$$ L_h(f) = \frac{h}{2} \left( f(x_0) + 2 \sum_{j=1}^{n-1} f(x_j) + f(x_n) \right), \quad (23) $$

where $x_j = a + j h$ with $h = (b-a)/n$. $L_h(f)$ is the approximation formula for the composite trapezoidal rule.

1. Richardson extrapolation can be used with the trapezoidal rule.

2. The value of $L_h(f)$ can be used to reduce the number of operations in the computation of $L_{h/2}(f)$.

Suppose that $L_h(f)$ is given by (23). Let $\tilde{n} = 2n$, $\tilde{h} = \frac{b-a}{\tilde{n}} = \frac{h}{2}$ and $\tilde{x}_j = a + j \tilde{h}$. Then, $\tilde{x}_{2j} = x_j$ for all $j$, leads to:

$$ 
\begin{aligned}
L_{h/2}(f) &= \frac{\tilde{h}}{2} \left( f(\tilde{x}_0) + 2 \sum_{j=1}^{\tilde{n}-1} f(\tilde{x}_j) + f(\tilde{x}_{\tilde{n}}) \right) = \frac{1}{2} \left( \frac{h}{2} \left( f(x_0) + 2 \sum_{j=1}^{n-1} f(x_j) + f(x_n) \right) + h \sum_{j=1}^{n} f(\tilde{x}_{2j-1}) \right) \\
&= \frac{1}{2} \left( L_h(f) + h \sum_{j=1}^{n} f(\tilde{x}_{2j-1}) \right),
\end{aligned} 
\quad (24) $$

### H. Adaptive Quadrature Methods

Let $f : [a, b] \to \mathbb{R}$ be a sufficiently continuously differentiable function. The goal is to approximate the integral

$$ \int_a^b f(x) dx $$

to a chosen accuracy.

In the composite methods, the step (the distance between $x_i$ partitioning the interval $[a, b]$) is constant. To get a good approximation of the integral, it would be advantageous to choose a smaller step where the function $f$ varies more rapidly. This is the idea motivating the adaptive quadrature methods.
There are many adaptive quadrature methods. The adaptive quadrature method considered in this part is based on the composite Simpson rule.

Let $a = x_0 < x_1 < \dots < x_n = b$ be a partition of $[a, b]$. The $x_i$ may not be equally spaced. Computing two approximations of

$$I_i = \int_{x_{i-1}}^{x_i} f(x) dx$$

using the composite Simpson rule. With $m = 1$, $a = x_{i-1}$, $b = x_i$ and $h = h_i = (x_i - x_{i-1})/2$ in Theorem 3, yields $I_i = S_i + R_i$, where

$$S_i = \frac{h_i}{3} (f(x_{i-1}) + 4f(x_{i-1} + h_i) + f(x_i)), \quad R_i = -\frac{f^{(4)}(\eta_i) (x_i - x_{i-1})}{180} h_i^4 = -\frac{f^{(4)}(\eta_i)}{90} h_i^5 \quad (25)$$

for $\eta_i$ between $x_{i-1}$ and $x_i$. Taking $m = 2$, $a = x_{i-1}$, $b = x_i$ and $h = h_i/2 = (x_i - x_{i-1})/4$ in Theorem 3, gives $I_i = \tilde{S}_i + \tilde{R}_i$, where

$$\tilde{S}_i = \frac{h_i}{6} \left( f(x_{i-1}) + 4f \left( x_{i-1} + \frac{h_i}{2} \right) + 2f(x_{i-1} + h_i) + 4f \left( x_{i-1} + \frac{3h_i}{2} \right) + f(x_i) \right), \\ \tilde{R}_i = -\frac{f^{(4)}(\mu_i) (x_i - x_{i-1})}{180} \left( \frac{h_i}{2} \right)^4 = -\frac{f^{(4)}(\mu_i)}{90 \times 16} h_i^5, \quad (26)$$

for $h_i = (x_i - x_{i-1})/2$ and $\mu_i$ between $x_{i-1}$ and $x_i$.

**Note:**

A very accurate code called CADRE (Cautious Adaptive Romberg Extrapolation) has been developed by Carl De Boor [2, 22] and integrated in IMSL [23] software package. This code is exploited in Section III.

## I. Filon method

Filon quadrature [6] is widely used in physics and engineering for computation of Fourier-type integrals. Applications include evaluation of highly oscillatory Sommerfeld integrals [5] encountered in electromagnetic and seismic problems in layered media and numerical solution of fluid mechanical problems, as well as various problems of neutron scattering and quantum mechanics.

The method is applied to oscillatory definite integrals of the form:

$$C = \int_{a}^{b} f(x) \cos(m\pi x) dx, \quad S = \int_{a}^{b} f(x) \sin(m\pi x) dx, \quad (27)$$

where $f(x)$ is a relatively slowly-varying function and $m$ is an integer. The domain $[a, b]$ is divided into $2N$ intervals of length $h$, where interpolation by parabolas is used. Since each interval is converted into a Fourier integral of quadratic polynomials, it is possible to obtain a closed-form via integration by parts.

The integration formulas are given as (taking $k = m\pi$):

$$C = \int_{a}^{b} f(x) \cos(kx) dx \approx h(\alpha[f(b)\sin(kb) - f(a)\sin(ka)] + \beta C_{2n} + \gamma C_{2n-1}) \\ S = \int_{a}^{b} f(x) \sin(kx) dx \approx h(\alpha[f(b)\cos(kb) - f(a)\cos(ka)] + \beta S_{2n} + \gamma S_{2n-1}) \quad (28)$$

where

$$\begin{aligned} \alpha &= \frac{(\theta^2 + \theta \sin(\theta) \cos(\theta) - 2 \sin^2(\theta))}{\theta^3} \\ \beta &= \frac{2[\theta(1 + \cos^2(\theta)) - 2 \sin(\theta) \cos(\theta)]}{\theta^3} \\ \gamma &= \frac{4(\sin(\theta) - \theta \cos(\theta))}{\theta^3}, \quad \theta = kh \end{aligned} \quad (29)$$

and:

$$\begin{aligned} C_{2n} &= \frac{1}{2} f(a) \cos(ka) + f(a + 2h) \cos(k(a + 2h)) + f(a + 4h) \cos(k(a + 4h)) + \dots + \frac{1}{2} f(b) \cos(kb) \\ C_{2n-1} &= f(a + h) \cos(k(a + h)) + f(a + 3h) \cos(k(a + 3h)) + \dots + f(b - h) \cos(k(b - h)), \\ S_{2n} &= \frac{1}{2} f(a) \cos(ka) + f(a + 2h) \sin(k(a + 2h)) + f(a + 4h) \cos(k(a + 4h)) + \dots + \frac{1}{2} f(b) \sin(kb) \\ S_{2n-1} &= f(a + h) \sin(k(a + h)) + f(a + 3h) \cos(k(a + 3h)) + \dots + f(b - h) \sin(k(b - h)) \end{aligned} \quad (30)$$
The error of Filon method [6] is given by $(b-a)MH(\theta)h^3 + O(h^4)$ where:

$$M = \max_{a \le x \le b} |f'''(x)|, \quad H(\theta) = \left| \frac{\sin \theta}{3\theta^2} + \frac{\cos \theta}{\theta^3} - \frac{\sin \theta}{\theta^4} \right| \quad (31)$$

## J. Monte-Carlo quadrature

In contrast to several deterministic methods described previously, a probabilistic method such as Monte-Carlo (MC) methods can be used also in quadrature since it has several advantages such as:

* Solving Dimension Curse [2]

In Quadrature problems, error is extremely dependent on system dimension (cf. Table 2). As shown in Table 2,

<table><caption>**Table 2:** Integration algorithms and associated errors in 1 and D dimensions with n subdivisions (in deterministic quadrature) or steps (in Monte-Carlo). MC error is ~ 1/√n. In Gauss method, m is the Legendre polynomial degree [2, 14] used for integration.</caption><thead><tr><th>Quadrature</th><th>Error in 1 dim.</th><th>Error in D dim.</th></tr></thead><tbody><tr><td>Trapezoidal</td><td>1/n<sup>2</sup></td><td>1/n<sup>2/D</sup></td></tr><tr><td>1/3 and 3/8 Simpson</td><td>1/n<sup>4</sup></td><td>1/n<sup>4/D</sup></td></tr><tr><td>Gauss</td><td>1/n<sup>2m-1</sup></td><td>1/n<sup>2m-1/D</sup></td></tr><tr><td>Monte-Carlo</td><td>1/√n</td><td>1/√n</td></tr></tbody></table>

when system dimension D grows, error increases since one has to use more meshing points as D increases whereas in MC, it is independent of D.

* Limited computational resources:

Since error is given by the number of steps, it suffices to let a computer with limited resources run for a time long enough to achieve the required steps. Thus runtime length compensates lack of resources.

* Versatility: MC methods can be used in many branches of Physics and Engineering such as Mechanics, Electro-magnetism, Statistical Physics, Quantum Mechanics, Atomic and Nuclear Physics, Field Theory, High Energy Physics... In numerical analysis, MC methods are used beside Quadrature, in Linear Algebra, root finding of vector functions, Optimization, Ordinary and Partial Differential Equations...

* MC methods can be improved in terms of accuracy and speed with Variance reduction techniques [13, 24]

In MC quadrature the sampling points $x_i$ are random variables that are picked from some probability distribution and the function $f(x)$ is thus sampled in a random fashion instead of a deterministic one [14]. The naive way to evaluate an integral in a random fashion is to express it as:

$$\int_a^b dxf(x) \approx \frac{1}{n} \sum_{i=1}^n f(x_i) \quad (32)$$

However there is a major paradigm shift in MC, in the way quadrature is performed since random points $(x_i, y_i)$ are generated in order to evaluate the area under $f(x)$ as shown in Fig. 11.

In order to perform integration of $f(x)$ shown in Fig. 11 function values are bound over $[a, b]$ and random points $x_i, y_i$ are generated such that $x_i \sim a + (b-a)U(a, b)$ and $y_i \sim U(0, c)$ where $U(a, b)$ (resp. $U(0, c)$) is a uniform random generator [13, 24] over the interval $[a, b]$ (resp. $[0, c]$). Since $y_i$ is uniformly distributed over $[0, c]$, its position is checked with the condition that $y_i \le f(x_i)$.

Then for $i = 1...n$, a counter count stores the number of times, a point $(x_i, y_i) \in D$.

The probability of being in $D$ is the ratio of the surface below $f(x)$ i.e. the integral $\int_a^b dxf(x)$ and the rectangular domain $[a, b] \times [0, c]$ surface.

Moreover the probability of being in $D$ is the total number $n$ of being in $D$ over the total number of trials $n$. Thus $\frac{\text{count}}{n} = \frac{S}{(b-a)c}$ and the integral value is $S = \frac{\text{count}(b-a)c}{n}$.

Monte Carlo quadrature error is estimated after defining the statistical quantity $\sigma_n$ [13, 25, 26]:

$$\sigma_n^2 = \frac{1}{n} \sum_{i=1}^{n} f^2(x_i) - \left[ \frac{1}{n} \sum_{i=1}^{n} f(x_i) \right]^2. \quad (33)$$
<figure><img src="image_14.png" /><figcaption>**Fig.11:** Function and random meshing at $x_i$, with $i = 0, \dots, n$, where $x_0 = a$ and $x_n = b$. $f(x)$ is bound by $c$ over $[a, b]$. The surface of domain $\mathcal{D}$ is bound by $[a, b]$ and the function $f(x)$ over $[a, b]$ is the integral $\int_a^b dx f(x)$. Integration is done by throwing random points $(x_i, y_i)$ and counting the number of times they fall into $\mathcal{D}$.</figcaption></figure>

The quantity $\pm\sigma_n/\sqrt{n}$ is the statistical error.

This explains the discrepancy obtained in the evaluation of the rectangle solid angle: $\sigma_n$ pertaining to function given by eq. 11 was too large (cf. Fig. 7) and when the transformed function eq. 12 was used instead, perfect agreement was obtained between MC evaluation and the analytical result (cf. Fig. 8).

As a conclusion, in order to obtain tangible results from MC quadrature, $n$ should be large and $\sigma_n$ should be small i.e. variance should be reduced as explained next.

### Variance reduction with importance sampling

The mean value $\mu$ of a function $f(X)$, with $X$ a random variable is given by:

$$ \mu = \int_{-\infty}^{+\infty} f(x)p(x)dx, \qquad (34) $$

$p(x)$ is the probability density function of the random variable $X$. In MC, a random sequence $\{x_i\}, i = 1, 2, \dots, n$ is sampled from $p(x)$ to yield:

$$ \langle f \rangle_n = \frac{1}{n} \sum_{i=1}^{n} f(x_i). \qquad (35) $$

In the limit $n \to \infty$, $\langle f \rangle_n \to \mu$. Moreover, by Central Limit Theorem [13, 26], when $n \to \infty$, the probability density of $\langle f \rangle_n$ tends to a Gaussian (cf. eq. 37) with mean $\mu$ and variance $\sigma^2/n$, where

$$ \sigma^2 = \int_{-\infty}^{+\infty} [f(x) - \mu]^2 p(x)dx \qquad (36) $$

MC estimate of $\mu$ is given by $\langle f \rangle_n \pm \sigma/\sqrt{n}$. This implies that $\langle f \rangle_n$ is expected to lie within $\pm\sigma/\sqrt{n}$ around $\mu$ (for $n$ sufficiently large), with a probability $p_n$ given by:

$$ p_n = \frac{\sqrt{n}}{\sigma\sqrt{2\pi}} \int_{\mu-(\sigma/\sqrt{n})}^{\mu+(\sigma/\sqrt{n})} \exp\left[-\frac{n(x-\mu)^2}{2\sigma^2}\right] dx \approx \frac{1}{\sqrt{2\pi}} \int_{-1}^{+1} \exp\left[-\frac{x^2}{2}\right] dx = 0.6826895 \qquad (37) $$

In fact, $\sigma$ is not known; consequently it is approximated by its MC estimate $\sigma_n$ given by eq. 33.

$\pm\sigma_n/\sqrt{n}$ is the resulting statistical error. Notice that sample size $n$ must be large enough for the above estimates of error and mean to be acceptable. Another valuable test is to check whether sample mean has reached a Gaussian distribution.

The error is of order $1/\sqrt{n}$ with $n$ the number of MC runs. While an MC method allows to avoid dimension curse, it is extremely slow otherwise, since for a reasonable accuracy such as $10^{-6}$, the number of steps $n$ should be about $10^{12}$ when $\sigma_n$ is ignored.
The remedy is to resort to variance reduction techniques [24] that minimize variance without altering the mean. The basic principle behind these techniques is called importance sampling [13, 14] and the entire procedure is named "importance-based".

Importance sampling helps select from sample space, regions providing larger contributions, instead of those leading to small irrelevant values to the MC sum. Suppose {$x_i$}, $i = 1, 2, \dots, n$ are sampled from an importance density $g(x)$ instead of $p(x)$. In order to preserve the mean, we define a modified function $F(x)$ given by:

$$F(x) = f(x) \frac{p(x)}{g(x)}. \qquad (38)$$

When $F$ mean is evaluated over the importance density $g(x)$, it is equal to the mean of the original function $f(x)$ evaluated over density $p(x)$:

$$\mu(F) = \int_{-\infty}^{+\infty} F(x)g(x)dx = \int_{-\infty}^{+\infty} f(x)p(x)dx = \mu(f) \qquad (39)$$

Thus {$x_i$}, $i = 1, 2, \dots, n$ sampled from $g(x)$, allow to calculate $p(x_i)/g(x_i)$. The importance-based estimate of $\mu$ is given by:

$$\langle F \rangle_n = \frac{1}{n} \sum_{i=1}^{n} f(x_i) \frac{p(x_i)}{g(x_i)}. \qquad (40)$$

When $n \to \infty$, $\langle F \rangle_n \to \mu$.

It is now possible to calculate the statistical error associated with $\langle F \rangle_n$. The second moment [14] is a valuable measure, since mean is preserved under importance sampling [13].

The importance-based second moment [24] using $g(x)$ is given by:

$$\langle F^2 \rangle_g = \int_{-\infty}^{+\infty} F^2(x)g(x)dx = \int_{-\infty}^{+\infty} \left[ \frac{p(x)}{g(x)} \right] f^2(x)p(x)dx \qquad (41)$$

The ordinary second moment [24] using $p(x)$ is given by:

$$\langle f^2 \rangle_p = \int_{-\infty}^{+\infty} f^2(x)p(x)dx \qquad (42)$$

Thus choosing an adequate importance function $g(x)$ and ensuring that, on average, the ratio $[p(x)/g(x)] \ll 1$, it is possible to make $\langle F^2 \rangle_g \ll \langle f^2 \rangle_p$. Consequently, sampling from an importance density helps estimate the mean with a much better statistical accuracy for a given sample size, since the corresponding variance has been reduced.

### III. NUMERICAL EXAMPLES

#### A. Estimation of three special integrals

$$\int_{0}^{1} \exp(-x^2) dx, \int_{0}^{\pi/2} \sin^2(x) dx, \int_{0}^{\infty} \frac{\exp(-x)}{\ln(x)} dx \qquad (43)$$

Exact values: $\int_{0}^{1} \exp(-x^2) dx = \frac{\sqrt{\pi}}{2} \operatorname{erf}(1)$, $\int_{0}^{\pi/2} \sin^2 x dx = \frac{\pi}{4}$, $PV \int_{0}^{\infty} \frac{\exp(-x)}{\ln(x)} dx = -0.1544796$ where $PV$ is Principal Value (real counterpart of Cauchy Principal Value).

The first integral can be evaluated analytically using the fact the Error function $\operatorname{erf}(x) = \frac{2}{\sqrt{\pi}} \int_{0}^{x} \exp(-y^2) dy$ yielding $\int_{0}^{1} \exp(-x^2) dx = \frac{\sqrt{\pi}}{2} \operatorname{erf}(1) \sim 0.746824$.
<figure><img src="image_15.png" /><figcaption>Fig.12: Graphs of functions to integrate: $\exp(-x^2)$, $\sin^2(x)$ et $\frac{\exp(-x)}{\ln(x)}$. Last function is singular at $x = 1$ and must be evaluated by its principal value.</figcaption></figure>

The second integral can be written as $\sin^2 x = \frac{1}{2}[1 - \cos(2x)]$, yielding $\int_{0}^{\pi/2} [1 - \cos(2x)] dx = \frac{\pi}{4} \sim 0.785399$.

The last integral is considered improper (due to a singularity of the integrand $\frac{\exp(-x)}{\ln(x)}$ at $x = 1$). It is possible to estimate its PV (Principal Value) through avoiding $x = 1$ by $\epsilon$ surrounding the singularity:

$$ \lim_{\epsilon \to 0} \int_{0}^{1-\epsilon} + \int_{1+\epsilon}^{\infty} \frac{\exp(-x)}{\ln(x)} dx = -2.649109 + 2.494639 = -0.1544704, \quad (44) $$

with $\epsilon = 10^{-3}$. Validation of the result is performed by taking several $\epsilon$ values. With $\epsilon = 10^{-4}$, -0.1544219 is obtained and for $\epsilon = 10^{-5}$, -0.1544264 both in agreement with Gradstein-Ryzhik [8] value: -0.1544796.

## B. Additional examples

1. Evaluate the following integrals to 6 digits by trapeze method.

Trapeze method has error $\frac{1}{N^2}$ with $N$ uniform divisions of $[a, b]$; thus $N \sim 10^3$ in order to have $10^{-6}$ error.

Hence $\int_{0}^{1} \exp(-x^2) dx$ is to evaluate with $2^{(n-1)}$ points. It suffices $n = 11$ to get $N = 2^{(n-1)} \sim 10^3$.

Numerical value of integral: 0.746824, Exact value: 0.746824

The second integral yields: $\int_{0}^{\pi/2} \sin^2 x dx = \frac{\pi}{4}$ with $2^{(n-1)}$ points

Numerical value of integral: 0.785398, Exact value: 0.785398

2. Evaluate the following integrals to 6 digits by Simpson method.

Simpson method has $\frac{1}{N+1}$ error with $N$ uniform divisions of $[a, b]$. In order to get $10^{-6}$ error, $N \sim \sqrt[4]{10^6} \sim 32$.

0.746824 is obtained for $\int_{0}^{1} \exp(-x^2) dx$ whereas $\int_{0}^{\pi/2} \sin^2 x dx$ yields 0.785398.

3. Evaluate the following integrals to 6 digits by Gauss-Legendre method.

Order $n$ Gauss-Legendre method is of error $\frac{f^{(2n)}(\xi)}{2n!}$ with $\xi \in [a, b]$. An error $10^{-6}$ requires $2n! \sim 10^6$, i.e. $n \sim 5$.

Thus: $\int_{0}^{1} \exp(-x^2) dx$ requires 10 Gauss-Legendre points.

Numerical value of integral: 0.746824, Exact value: 0.746824

For the second integral: $\int_{0}^{\pi/2} \sin^2 x dx = \frac{\pi}{4}$ with 10 Gauss-Legendre points.

Numerical value of integral: 0.785398, Exact value: 0.785398
4. Evaluate the following integral to 16 digits:

$$ \int_{1}^{2} \frac{dx}{x} = \ln 2 = 0.69314\ 71805\ 59945\ 31 \qquad (45) $$

This integral can be evaluated with $\ln 2$ calculated in double precision or using CADRE (Cautious ADaptive Romberg Extrapolation) developed by Carl De Boor [2, 22, 27] or with multiple precision based on manipulating blocks of numbers like symbols as done by David H. Bailey ARPREC (ARbitrary PRECision) software at University of California, Davis and Lawrence Berkeley National Laboratory.

Using Simpson method ($\frac{1}{N^4}$ error when $N$ uniform divisions are performed over $[a, b]$) in double precision, $N \sim 10^4$ are needed to obtain $10^{-16}$ error. The result: 0.69314 71805 59944 8, with $4.10^{-16}$ error validate the precision required.

## APPENDIX A: NEWTON-COTES FORMULAS

Starting with local integration:

$$ \int_{a}^{a+mh} f(x)dx = \frac{1}{D}[f(a+mh) - f(a)] = \frac{h[(1+\Delta)^m - 1]}{\ln(1+\Delta)}f(a) \qquad (A1) $$

$\Delta$ is the finite difference operator whereas $\frac{1}{D}$ is the continuous integration operator [27]. Replacing $\Delta$ by $x$ and performing Taylor expansion around $x = 0$ yields:

$$ \frac{(1+x)^m - 1}{\ln(1+x)} \approx m + \frac{m^2 x}{2} + \frac{(2m^3 - 3m^2)x^2}{12} + \frac{(m^4 - 4m^3 + 4m^2)x^3}{24} + \frac{(6m^5 - 45m^4 + 110m^3 - 90m^2)x^4}{720} + \frac{(2m^6 - 24m^5 + 105m^4 - 200m^3 + 144m^2)x^5}{1440} + \frac{(12m^7 - 210m^6 + 1428m^5 - 4725m^4 + 7672m^3 - 5040m^2)x^6}{60480} \qquad (A2) $$

This is understood by performing the Taylor expansion of $\frac{(1+x)^m - 1}{\ln(1+x)}$ to order $m + 1$, substituting $m$ successively from 1 to 7 as displayed in Table 3:

<table><caption>**Table 3:** Taylor expansion of $(\frac{1+x)^m - 1}{\ln(1+x)}$ to order $m + 1$ obtained with MAXIMA.</caption><thead><tr><th>Order m</th><th>Expansion to order m + 1</th></tr></thead><tbody><tr><td>m = 1</td><td>1 + x/2 + x<sup>2</sup>/12 + ...</td></tr><tr><td>m = 2</td><td>2 + 2x + x<sup>2</sup>/3 + ...</td></tr><tr><td>m = 3</td><td>3 + 9x/2 + 9x<sup>2</sup>/4 + 3x<sup>3</sup>/8 + 3x<sup>4</sup>/80 + ...</td></tr><tr><td>m = 4</td><td>4 + 8x + 20x<sup>2</sup>/3 + 8x<sup>3</sup>/3 + 14x<sup>4</sup>/45 + ...</td></tr><tr><td>m = 5</td><td>5 + 25x/2 + 175x<sup>2</sup>/12 + 75x<sup>3</sup>/8 + 425x<sup>4</sup>/144 + 95x<sup>5</sup>/288 + 275x<sup>6</sup>/12096 + ...</td></tr><tr><td>m = 6</td><td>6 + 18x + 27x<sup>2</sup> + 24x<sup>3</sup> + 123x<sup>4</sup>/10 + 33x<sup>5</sup>/10 + 41x<sup>6</sup>/140 + ...</td></tr><tr><td>m = 7</td><td>7 + 49x/2 + 539x<sup>2</sup>/12 + 1225x<sup>3</sup>/24 + 26117x<sup>4</sup>/720 + 2499x<sup>5</sup>/160 + 30919x<sup>6</sup>/8640 + 5257x<sup>7</sup>/17280 + 8183x<sup>8</sup>/518400 + ...</td></tr></tbody></table>

### 1. Open and Closed formulas

There are two types of Newton-Cotes formulas: closed type uses the function value at all points in $[a, b]$, whereas open type uses points in the interval $]a, b[$.

Closed Newton-Cotes formula of degree $n$ is stated as: $\int_a^b f(x) dx \approx \sum_{i=0}^n w_i f(x_i)$ where $x_i = i.h + x_0$, with $h$ the step equal to $(x_n - x_0)/m = (b - a)/n$.
The weights $w_i$ are obtained from Lagrange polynomials. They depend only on the $x_i$ and not on the function $f$. Let $L(x)$ be the Lagrange interpolation polynomial for the given data points $(x_0, f(x_0)), \dots, (x_n, f(x_n))$, then
$$ \int_a^b f(x) dx \approx \int_a^b L(x) dx = \int_a^b \left( \sum_{i=0}^n f(x_i) l_i(x) \right) dx $$

The open Newton-Cotes formula of degree $n \neq 0$ is $\int_a^b f(x) dx \approx \sum_{i=1}^{n-1} w_i f(x_i)$ with $w_i = \int_a^b l_i(x) dx$
The weights are found in a manner akin to the closed case.

## 2. Special case $m = 0$

For $m = 0$, Newton-Cotes approach is not valid. Starting from the finite step expansion: $f(x) = f(a) + (x-a)f'(\xi)$, the integration error is:

$$ \epsilon = \int_{a}^{b} [f(x) - P_{0}(x)] dx = \int_{a}^{b} [f(x) - f(a)] dx = \int_{a}^{b} (x-a)f'(\xi) dx = f'(\xi) \int_{a}^{b} (x-a) dx \\ = \frac{(b-a)^2}{2} f'(\xi) = \frac{h^2}{2} f'(\xi), \quad \xi \in [a, b] \qquad (A3) $$

## 3. Closed Newton-Cotes formulas: $m \neq 0$

This table lists some of the Newton-Cotes formulas of the closed type. For $0 \le i \le m$, with $m$ the degree, let $x_i = a + i \frac{b-a}{n} = a + ih$, and the notation $f_i$ be a shorthand for $f(x_i)$.

<table><caption>**Table 4:** Closed Newton-Cotes formulas for local errors.</caption><thead><tr><th>Degree n</th><th>Step size h</th><th>Common names</th><th>Formula</th><th>Error term</th></tr></thead><tbody><tr><td>1</td><td>h = (b - a)</td><td>Trapezoid rule</td><td><sup>h</sup>&frasl;<sub>2</sub>(f<sub>0</sub> + f<sub>1</sub>)</td><td>-<sup>(b-a)3</sup>&frasl;<sub>12</sub> f<sup>(2)</sup>(ξ)</td></tr><tr><td>2</td><td>h = <sup>(b-a)</sup>&frasl;<sub>2</sub></td><td>1/3 Simpson rule</td><td>(1/3) &middot; h(f<sub>0</sub> + 4f<sub>1</sub> + f<sub>2</sub>)</td><td>-<sup>(b-a)5</sup>&frasl;<sub>90.25</sub> f<sup>(4)</sup>(ξ)</td></tr><tr><td>3</td><td>h = <sup>(b-a)</sup>&frasl;<sub>3</sub></td><td>3/8 Simpson rule</td><td>(3/8) &middot; h(f<sub>0</sub> + 3f<sub>1</sub> + 3f<sub>2</sub> + f<sub>3</sub>)</td><td>-<sup>3(b-a)5</sup>&frasl;<sub>80.35</sub> f<sup>(4)</sup>(ξ)</td></tr><tr><td>4</td><td>h = <sup>(b-a)</sup>&frasl;<sub>4</sub></td><td>Boole rule</td><td><sup>2h</sup>&frasl;<sub>45</sub>(7f<sub>0</sub> + 32f<sub>1</sub> + 12f<sub>2</sub> + 32f<sub>3</sub> + 7f<sub>4</sub>)</td><td>-<sup>8(b-a)7</sup>&frasl;<sub>945.47</sub> f<sup>(6)</sup>(ξ)</td></tr></tbody></table>

The exponent of the segment size $(b-a)$ in the error term shows the rate at which the approximation error decreases. The degree of the derivative of $f$ in the error term gives the degree up to which polynomials can be integrated exactly with this rule. Note that the derivative of $f$ in the error term increases by 2 for every other rule and $\xi \in [a, b]$.

## 4. Open Newton-Cotes formulas: $m \neq 0$

This table lists some of the Newton-Cotes formulas of the open type. Again, $f_i$ is a shorthand for $f(x_i)$, with $x_i = a + i \left(\frac{b-a}{n}\right)$, and the degree.

<table><caption>**Table 5:** Open Newton-Cotes formulas for local errors.</caption><thead><tr><th>Common name</th><th>step size h</th><th>Formula</th><th>Error term</th><th>Degree</th></tr></thead><tbody><tr><td>Rectangle/Midpoint rule</td><td><sup>b-a</sup>&frasl;<sub>2</sub></td><td>2hf<sub>1</sub></td><td><sup>1</sup>&frasl;<sub>3</sub>h<sup>3</sup>f<sup>(2)</sup>(ξ)</td><td>2</td></tr><tr><td>Trapezoid method</td><td><sup>b-a</sup>&frasl;<sub>3</sub></td><td><sup>3</sup>&frasl;<sub>2</sub>h(f<sub>1</sub> + f<sub>2</sub>)</td><td><sup>1</sup>&frasl;<sub>4</sub>h<sup>3</sup>f<sup>(2)</sup>(ξ)</td><td>3</td></tr><tr><td>Boole rule</td><td><sup>b-a</sup>&frasl;<sub>4</sub></td><td><sup>4</sup>&frasl;<sub>3</sub>h(2f<sub>1</sub> - f<sub>2</sub> + 2f<sub>3</sub>)</td><td><sup>28</sup>&frasl;<sub>90</sub>h<sup>5</sup>f<sup>(4)</sup>(ξ)</td><td>4</td></tr><tr><td>Milne rule</td><td><sup>b-a</sup>&frasl;<sub>5</sub></td><td><sup>5</sup>&frasl;<sub>24</sub>h(11f<sub>1</sub> + f<sub>2</sub> + f<sub>3</sub> + 11f<sub>4</sub>)</td><td><sup>95</sup>&frasl;<sub>144</sub>h<sup>5</sup>f<sup>(4)</sup>(ξ)</td><td>5</td></tr></tbody></table>

## 5. Instability at high degree

A Newton-Cotes formula of any degree $n$ can be constructed. However, for large $n$ a Newton-Cotes rule may suffer from Runge phenomenon where the error grows exponentially with $m$. Methods such as Gaussian quadrature with
unequally spaced points are stable and more accurate than Newton-Cotes. When the integrand is given at a fixed uniform grid, these methods cannot be used, then Runge phenomenon can be avoided via the use of a composite rule.

APPENDIX B: QUADRATURE IN THE COMPLEX PLANE

While this work describes quadrature in a real 1D or higher dimensional domains, analytic continuation allows to perform quadrature by complex integration [28] leading to the sought result in a much easier way. Consider the integral:

$$I_n = \int_C (z - z_0)^n dz, \qquad (B1)$$

where C is a closed contour surrounding $z = z_0$ and n is an integer. Since $(z-z_0)^n$ is analytic, except possibly at $z_0$, we may deform the contour into a circle centered at $z_0$ without affecting $I_n$. Then we may write $z-z_0 = Re^{i\theta}$, $dz = iRe^{i\theta}d\theta$ to write:

$$I_n = iR^{n+1} \int_0^{2\pi} e^{i(n+1)\theta} d\theta = \begin{cases} 0, & n \neq -1 \\ 2i\pi, & n = -1 \end{cases} \qquad (B2)$$

Now suppose we have a function, $f(z)$, which is analytic in a region except at the point $z_0$ in the region. Then we can make the Laurent expansion [28] about $z_0$:

$$f(z) = \sum_{n=-\infty}^{\infty} A_n(z - z_0)^n. \qquad (B3)$$

Take a contour C around $z_0$:

$$\frac{1}{2\pi i} \int_C f(z) dz = \frac{1}{2\pi i} \int_C \sum_{n=-\infty}^{\infty} A_n(z - z_0)^n dz \qquad (B4)$$

$$= A_{-1}. \qquad (B5)$$

Thus the coefficient of $1/(z-z_0)$ in the Laurent series is given by $A_{-1} = \frac{1}{2\pi i} \int_C f(z)dz$. This coefficient is called the residue of $f(z)$ at $z_0$. Notice that the residue is zero if $f(z)$ is analytic at $z_0$, or if the coefficient $A_{-1}$ is zero (even if $z_0$ is a pole or isolated essential singularity).

An important and useful notion is the residue theorem [28]. Consider contour C in a region where $f(z)$ is analytic except at isolated singularities (poles or essential singularities) as in Fig. 13.

<figure><img src="image_16.png" /><figcaption>Fig.13:</figcaption></figure>

Contours to illustrate residue theorem. Singularities are at a, b, and c.

We want to determine $\int_C f(z) dz$. We can write this in terms of the sum of the integrals around each singularity:

$$\begin{aligned} \int_C f(z) dz &= \int_{C_a} f(z) dz + \int_{C_b} f(z) dz + \dots \\ &= 2i\pi(a_{-1} + b_{-1} + c_{-1} + \dots) \\ &= 2i\pi \sum_{a,b,c,\dots} R_e, \end{aligned} \qquad (B6)$$
where $a_{-1}, b_{-1}, c_{-1}, \dots$ are the residues at $a, b, c, \dots$, respectively, and $\sum R_e$ is the sum of the residues of $f(z)$ interior to the contour $C$.

The computation of the residues is thus often an important part of evaluating integrals. At a simple pole, the Laurent series is

$$f(z) = \frac{A_{-1}}{z - z_0} + \sum_{n=0}^{\infty} A_n (z - z_0)^n, \quad A_{-1} = \lim_{z \to z_0} [(z - z_0) f(z)]$$

For a pole of order $m$,

$$f(z) = \frac{A_{-m}}{(z - z_0)^m} + \frac{A_{-m+1}}{(z - z_0)^{m-1}} + \dots + \frac{A_{-1}}{z - z_0} + \sum_{n=0}^{\infty} A_n (z - z_0)^n.$$

If we multiply both sides by $(z - z_0)^m$, we have:

$$(z - z_0)^m f(z) = A_{-m} + A_{-m+1}(z - z_0) + \dots + A_{-1}(z - z_0)^{m-1} + \sum_{n=0}^{\infty} A_n (z - z_0)^{n+m}.$$

Now differentiate $m-1$ times and evaluate at $z = z_0$:

$$\left. \frac{d^{m-1} [(z - z_0)^m f(z)]}{dz^{m-1}} \right|_{z=z_0} = (m-1)! A_{-1},$$

and hence,

$$A_{-1} = \frac{1}{(m-1)!} \left. \frac{d^{m-1} [(z - z_0)^m f(z)]}{dz^{m-1}} \right|_{z=z_0}.$$

Sometimes it is easier to just carry out the Laurent expansion sufficiently to find $A_{-1}$ directly or consider a fast method for evaluating residues described as follows:

For a simple pole $z_0$, the residue [28] is $\lim_{z \to z_0} [(z - z_0) f(z)]$.

For an order $m > 1$ pole, supposing the integrand $f(z) = N(z)/D(z)$, then the residue [28] is $N(z_0)/D'(z_0)$, provided $D'(z_0) \neq 0$ without differentiating $f(z)$, $m-1$ times or performing a Laurent expansion around the $m > 1$ order pole.

<figure><img src="image_17.png" /><figcaption>**Fig.14:** Contour for complex integration of $\oint_C dz(z^{p-1})/(1+z)$ displaying a branch cut about the origin and a simple pole $z = -1$.</figcaption></figure>

Suppose $f(z)$ has a simple pole at $z = z_0 = x_0 + i0$ on the real axis. One may define an integral along the real axis through this pole according to: $PV \int_{\alpha}^{\beta} f(z) dz \equiv \lim_{\epsilon \to 0} \left[ \int_{\alpha}^{x_0 - \epsilon} f(x) dx + \int_{x_0 + \epsilon}^{\beta} f(x) dx \right]$, where $\alpha < x_0 < \beta$. This is Cauchy Principal Value Integral.

Consider evaluating the following integral [28]:

$$\int_{0}^{\infty} dx \frac{x^{p-1}}{1+x}, \quad 0 < p < 1 \qquad (B7)$$
Using analytic continuation, this is rewritten with $\oint_C dz(z^{p-1})/(1+z)$ where C is the contour shown in Fig. 14. Point $z=0$ is a branch point and must be avoided for functional definition unicity [28]. Additionally, there is a simple pole $z=-1$ with residue given by $\lim_{z \to -1}(z+1)(z^{p-1})/(1+z) = e^{(p-1)i\pi}$.

Performing contour integration along AB, BDEFG, GH, HJA path yields:

$$ \int_{\epsilon}^{R} dx \frac{x^{p-1}}{1+x} + \int_{0}^{2\pi} d\theta \frac{(Re^{i\theta})^{p-1} iRe^{i\theta}}{1+Re^{i\theta}} + \int_{R}^{\epsilon} dx \frac{(xe^{2i\pi})^{p-1}}{1+xe^{2i\pi}} + \int_{2\pi}^{0} d\theta \frac{(\epsilon e^{i\theta})^{p-1} i\epsilon e^{i\theta}}{1+\epsilon e^{i\theta}} = 2i\pi e^{(p-1)i\pi} \quad (B8) $$

Taking the limits $\epsilon \to 0, R \to \infty$ make the second and fourth integral vanish. Consequently, the final result is:

$$ \int_{0}^{\infty} dx \frac{x^{p-1}}{1+x} = \frac{p}{\sin p\pi} \quad (B9) $$

[1] L. Hodgkin, *A History of Mathematics From Mesopotamia to Modernity*, Oxford University Press Inc., New York (2005)

[2] J. R. Rice, *Numerical Methods, Software and Analysis*, Second Edition, Academic Press, New-York (1993).

[3] T. J. Osler, Vieta-like products of nested radicals with Fibonacci and Lucas numbers, The Fibonacci Quarterly, 52, 202 (2007) and The united Vieta and Wallis product for $\pi$, American Math. Monthly 106, 774 (1999).

[4] A. B. Migdal, *Qualitative methods in Quantum theory*, W. A. Benjamin, Inc. Advanced Book Program, Reading Mass. (1977).

[5] M. Abramowitz and I.S Stegun, *Handbook of Mathematical Tables*, Dover, New-York (1960). M. Abramowitz, *On the practical evaluation of integrals*, Journal of the Society for Industrial and Applied Mathematics (SIAM), Volume 2, (1954)

[6] P. J. Davis and P. Rabinowitz, *Methods of numerical integration*, Harcourt Brace Jovanovich, Publishers and Academic Press, Inc. New-York (1984)

[7] M. Tinkham, *Introduction to Superconductivity* (Dover Publications, New York, 2nd ed., 2004).

[8] I. S. Gradstein and I. M. Ryzhik, *Table of Integrals, Series and Products*, Fifth Edition, edited by Allan Jeffrey (Academic Press) (1994).

[9] J. W. Wilkins, *Solid State Physics Lectures*, Cornell University, Ithaca N.Y. (1982).

[10] D. Zwillinger, *CRC Standard Mathematical Tables and Formulae*, Chapman and Hall CRC, Boca Raton (2003).

[11] J. Casey, *A treatise on Spherical Trigonometry and its applications to Geodesy and Astronomy*, Hodges, Figgis and Co. Dublin (1889).

[12] F. Eriksson, *On the measure of solid angles*, Math. Mag. 63 (3): 184 (1990). A. Van Oosterom and J. Strakee, *The solid angle of a plane triangle IEEE Trans. on BioMed. Eng. BME-30*, 125 (1983)

[13] D. P. Landau and K. Binder, *A Guide to Monte-Carlo Simulations in Statistical Physics*, 2nd Edition, Cambridge (2000).

[14] W. H. Press, W. T. Vetterling, S. A. Teukolsky and B. P. Flannery, *Numerical Recipes in C: The Art of Scientific Computing*, Second Edition, Cambridge University Press (New-York, 1992).

[15] F. B. Daniels, *A Planimetric Method of Harmonic Analysis*, J. App. Phys. 23, 369 (1952)

[16] O. Mader, Elektrotech. Z. 36, 847 (1909).

[17] R. W. Gatterdam, *The planimeter as an example of Green theorem*, The American Mathematical Monthly, 88 (9): 701-704 (1981)

[18] G. Arfken, *Mathematical Methods for Physicists*, Academic Press, San Diego, (1985).

[19] Kelvin and Tait, *Natural Philosophy, Proc Roy. Soc, xxvii*, p.371 (1878), Henrici, Phil. Mag., xxxviii, p. 110, (1894), Sharp, Phil. Mag., xxxviii, p. 121, (1894), Yule, Phil. Mag., xxxix, p. 367 (1895), The Electrician, March 22 (1895), A. A. Michelson and S. W. Stratton, *A New Harmonic Analyzer*, American Journal of Science 25, 1-13 (1898)

[20] W. Hammack, S. Kranz and B. Carpenter, *A Visual Tour of a Nineteenth Century Machine that Performs Fourier Analysis Articulate Noise Books* (2014)

[21] B. Dionne, *MAT3380 and graduate courses in Numerical Analysis*, Department of Mathematics and Statistics, University of Ottawa (2023).

[22] S. D. Conte, C. de Boor, *Elementary Numerical Analysis*, McGraw-Hill, New York, (1980).

[23] IMSL Library, 8th edition. International Mathematical and Statistical Libraries, Inc., Houston, TX, (1980).

[24] J. E. Gentle, *Random number generation and Monte-Carlo methods*, Second Edition, Springer Science Business Media, Inc., New-York (2003).

[25] W. Feller, *An Introduction to Probability Theory and its Applications*, Vols. I and II, John Wiley, New York (1968).

[26] A. Papoulis, *Probability Theory, Random Variables and Stochastic Processes*, McGraw Hill (1965).

[27] J. Stoer R. and Bulirsch, *Introduction to Numerical Analysis*, Second Edition, Springer-Verlag, (New York, 1992).

[28] M. Ablowitz and A. Fokas, *Complex variables: Introduction and Applications*, 2nd edition, p. 287, Cambridge University Press, New-York (2003).