Yet Another

Introduction to Numerical Methods

    version 13.05

D.V. Fedorov

\begin{matrix}c \\ ◯ \end{matrix} 2013 Dmitri V. Fedorov

Permission is granted to copy and redistribute this work under the terms of either the GNU General Public License^{1}, version 3 or later, as published by the Free Software Foundation, or the Creative Commons Attribution Share Alike License^{2}, version 3 or later, as published by the Creative Commons corporation.

This work is distributed in the hope that it will be useful, but without any warranty. No responsibility is assumed by the author and the publisher for any damage from any use of any methods, instructions or ideas contained in the material herein.

1http://en.wikipedia.org/wiki/GPL

2http://en.wikipedia.org/wiki/CC-BY-SA

Preface

This book evolved from lecture notes developed over several years of teaching numerical methods at the University of Aarhus. The book contains short descriptions of some of the most common numerical methods together with illustrational implementations of the discussed algorithms mostly in the C programming language. The implementations are solely for educational purposes, they are not thoroughly tested and are not tuned in any way other than to fit nicely on a single page of the book; some simple optimizations were often sacrificed to conciseness. For serious calculations one should use the GNU Scientific Library.

The content of the book is free as in freedom. You are permitted to copy and redistribute the book in original or modified form either gratis or for a fee. However, you must attribute the original author(s) and pass the same freedom to all recipients of your copies. See the GPL and CC-BY-SA licenses for more details.

2013 Dmitri Fedorov iv

Contents

1 Interpolation 1 1.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Polynomial interpolation . . . . . . . . . . . . . . . . . . . . . 1 1.3 Spline interpolation . . . . . . . . . . . . . . . . . . . . . . . . 2

1.3.1 Linear interpolation . . . . . . . . . . . . . . . . . . . 3 1.3.2 Quadratic spline . . . . . . . . . . . . . . . . . . . . . 4 1.3.3 Cubic spline . . . . . . . . . . . . . . . . . . . . . . . . 4 1.3.4 Akima sub-spline interpolation . . . . . . . . . . . . . 7

1.4 Other forms of interpolation . . . . . . . . . . . . . . . . . . . 9 1.5 Multivariate interpolation . . . . . . . . . . . . . . . . . . . . 9

1.5.1 Nearest-neighbor interpolation . . . . . . . . . . . . . 9 1.5.2 Piecewise-linear interpolation . . . . . . . . . . . . . . 9 1.5.3 Bi-linear interpolation . . . . . . . . . . . . . . . . . . 9

2 Numerical integration 13 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 2.2 Rectangle and trapezium rules . . . . . . . . . . . . . . . . . 13 2.3 Quadratures with regularly spaced abscissas . . . . . . . . . . 15

2.3.1 Classical quadratures . . . . . . . . . . . . . . . . . . . 16 2.4 Quadratures with optimized abscissas . . . . . . . . . . . . . 17

2.4.1 Gauss quadratures . . . . . . . . . . . . . . . . . . . . 18 2.4.2 Gauss-Kronrod quadratures . . . . . . . . . . . . . . . 21

2.5 Adaptive quadratures . . . . . . . . . . . . . . . . . . . . . . 22 2.6 Variable transformation quadratures . . . . . . . . . . . . . . 23 2.7 Infinite intervals . . . . . . . . . . . . . . . . . . . . . . . . . 25

3 Monte Carlo integration 27 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 3.2 Plain Monte Carlo sampling . . . . . . . . . . . . . . . . . . . 28 3.3 Importance sampling . . . . . . . . . . . . . . . . . . . . . . . 28 3.4 Stratified sampling . . . . . . . . . . . . . . . . . . . . . . . . 29 3.5 Quasi-random (low-discrepancy) sampling . . . . . . . . . . . 31

3.5.1 Van der Corput and Halton sequences . . . . . . . . . 31

3.5.2 Lattice rules . . . . . . . . . . . . . . . . . . . . . . . 32 3.6 Implementations . . . . . . . . . . . . . . . . . . . . . . . . . 33

4 Ordinary differential equations 35 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 4.2 Error estimate . . . . . . . . . . . . . . . . . . . . . . . . . . 36 4.3 Runge-Kutta methods . . . . . . . . . . . . . . . . . . . . . . 37

4.3.1 Embeded methods with error estimates . . . . . . . . 39 4.4 Multistep methods . . . . . . . . . . . . . . . . . . . . . . . . 41

4.4.1 Two-step method . . . . . . . . . . . . . . . . . . . . . 41 4.4.2 Two-step method with extra evaluation . . . . . . . . 42

4.5 Predictor-corrector methods . . . . . . . . . . . . . . . . . . . 42

4.5.1 Two-step method with correction . . . . . . . . . . . . 43 4.6 Adaptive step-size control . . . . . . . . . . . . . . . . . . . . 43

5 Systems of linear equations 47 5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 5.2 Triangular systems . . . . . . . . . . . . . . . . . . . . . . . . 48 5.3 Reduction to triangular form . . . . . . . . . . . . . . . . . . 48

5.3.1 QR-decomposition . . . . . . . . . . . . . . . . . . . . 49 5.3.2 LU-decomposition . . . . . . . . . . . . . . . . . . . . 54

5.4 Determinant of a matrix . . . . . . . . . . . . . . . . . . . . . 55 5.5 Matrix inverse . . . . . . . . . . . . . . . . . . . . . . . . . . 55

6 Ordinary least squares 57 6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 6.2 Linear least-squares problem . . . . . . . . . . . . . . . . . . 57 6.3 Solution via QR-decomposition . . . . . . . . . . . . . . . . . 57 6.4 Ordinary least-squares curve fitting . . . . . . . . . . . . . . . 58

6.4.1 Variances and correlations of fitting parameters …. 59

7 Nonlinear equations 61 7.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 7.2 Newton’s method . . . . . . . . . . . . . . . . . . . . . . . . . 61 7.3 Broyden’s quasi-Newton method . . . . . . . . . . . . . . . . 63

8 Minimization and optimization 67 8.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 8.2 Local minimization . . . . . . . . . . . . . . . . . . . . . . . . 67

8.2.1 Newton’s methods . . . . . . . . . . . . . . . . . . . . 67 8.2.2 Quasi-Newton methods . . . . . . . . . . . . . . . . . 68 8.2.3 Downhill simplex method . . . . . . . . . . . . . . . . 69

8.3 Global optimization . . . . . . . . . . . . . . . . . . . . . . . 70 8.3.1 Simulated annealing . . . . . . . . . . . . . . . . . . . 71 8.3.2 Quantum annealing . . . . . . . . . . . . . . . . . . . 72 8.3.3 Evolutionary algorithms . . . . . . . . . . . . . . . . . 73

8.4 Implementation in C . . . . . . . . . . . . . . . . . . . . . . . 74

9 Eigenvalues and eigenvectors 77 9.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 9.2 Similarity transformations . . . . . . . . . . . . . . . . . . . . 77

9.2.1 Jacobi eigenvalue algorithm . . . . . . . . . . . . . . . 78 9.3 Eigenvalues of updated matrix . . . . . . . . . . . . . . . . . 80

9.3.1 Rank-1 update . . . . . . . . . . . . . . . . . . . . . . 81 9.3.2 Symmetric row/column update . . . . . . . . . . . . . 82 9.3.3 Symmetric rank-2 update . . . . . . . . . . . . . . . . 83

10 Power iteration methods and Krylov subspaces 85 10.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 10.2 Krylov subspaces . . . . . . . . . . . . . . . . . . . . . . . . . 86 10.3 Arnoldi iteration . . . . . . . . . . . . . . . . . . . . . . . . . 86 10.4 Lanczos iteration . . . . . . . . . . . . . . . . . . . . . . . . . 87 10.5 Generalised minimum residual (GMRES) . . . . . . . . . . . 88

11 Fast Fourier transform 89

11.1 Discrete Fourier Transform . . . . . . . . . . . . . . . . . . . 89

11.1.1 Applications . . . . . . . . . . . . . . . . . . . . . . . 90 11.2 Cooley-Tukey algorithm . . . . . . . . . . . . . . . . . . . . . 91 11.3 Multidimensional DFT . . . . . . . . . . . . . . . . . . . . . . 91 11.4 C implementation . . . . . . . . . . . . . . . . . . . . . . . . . 91 viii CONTENTS

Chapter 1

Interpolation

1.1 Introduction

In practice one often meets a situation where the function of interest, f(x), is only represented by a discrete set of tabulated points,

{x_{i}, y_{i} = f(x_{i}) | i = 1 … n},

obtained, for example, by sampling, experimentation, or extensive numerical calculations.

Interpolation means constructing a (smooth) function, called interpolating function or interpolant, which passes exactly through the given points and hopefully approximates the tabulated function in between the tabulated points. Interpolation is a specific case of curve fitting in which the fitting function must go exactly through the data points.

The interpolating function can be used for different practical needs like estimating the tabulated function between the tabulated points and estimating the derivatives and integrals involving the tabulated function.

1.2 Polynomial interpolation

Polynomial interpolation uses a polynomial as the interpolating function. Given a table of n points, {x_{i}, y_{i}}, where no two x_{i} are the same, one can construct a polynomial P (x) of the order n−1 which passes exactly through the points: P (x_{i}) = y_{i}. This polynomial can be intuitively written in the Lagrange form,

P (x) = ∑_{i=1}^{n} y_{i} ∏_{k≠i}^{n} \frac{x − x_{k}}{x_{i} − x_{k}} .  (1.1)

The interpolating polynomial always exists and is unique.

Higher order interpolating polynomials are susceptible to the Runge’s phenomenon: erratic oscillations close to the end-points of the interval (see Figure 1.1).

Table 1.1: Polynomial interpolation in C

double p o l i n t e r p ( int n , double ∗x , double ∗y , double z ) {
double s =0;
for ( int i =0; i<n ; i++) {
double p=1;
for ( int k=0;k<n ; k++) i f ( k!= i ) p∗=(z−x [ k ] ) / ( x [ i ]−x [ k ] ) ;
s+=y [ i ] ∗ p ; }
return s ; }

Figure 1.1: Lagrange interpolating polynomial, solid line, showing the Runge’s phenomenon: large oscillations at the edges. For comparison the dashed line shows a cubic spline.

1.3 Spline interpolation

Spline interpolation uses a piecewise polynomial, S(x), called spline, as the interpolating function,

S(x) = S_{i}(x) , if x ∈ [x_{i}, x_{i+1}] , i = 1, …, n − 1  (1.2)

where S_{i}(x) is a polynomial of a given order k.

The spline of the order k ≥ 1 can be made continuous at the tabulated points,

S_{i}(x_{i}) = y_{i} , S_{i}(x_{i+1}) = y_{i+1} , i = 1, …, n − 1 ,  (1.3)

together with its k − 1 derivatives,

S′_{i}(x_{i+1}) = S′_{i+1}(x_{i+1}) , ; S′′_{i}(x_{i+1}) = S′′_{i+1}(x_{i+1}) , ; ⋮ ; S_{i}^{(k−1)}(x_{i+1}) = S_{i+1}^{(k−1)}(x_{i+1}) . | i = 1, …, n − 2  (1.4)

Continuity conditions (1.3) and (1.4) make kn + n − 2k linear equations for the (n − 1)(k + 1) = kn + n − k − 1 coefficients of n − 1 polynomials (1.2) of the order k. The missing k − 1 conditions can be chosen (reasonably) arbitrarily.

The most popular is the cubic spline, where the polynomials S_{i}(x) are of third order. The cubic spline is a continuous function together with its first and second derivatives. The cubic spline has a nice feature that it (sort of) minimizes the total curvature of the interpolating function. This makes the cubic splines look good.

Quadratic splines—continuous with the first derivative—are not nearly as good as cubic splines in most respects. In particular they might oscillate unpleasantly when a quick change in the tabulated function is followed by a period where the function is nearly a constant. Cubic splines are somewhat less susceptible to such oscillations.

Linear spline is simply a polygon drawn through the tabulated points.

1.3.1 Linear interpolation

Linear interpolation is a spline with linear polynomials. The continuity conditions (1.3) can be satisfied by choosing the spline in the (intuitive) form

S_{i}(x) = y_{i} + p_{i}(x − x_{i}) ,  (1.5)

where

p_{i} = \frac{Δy_{i}}{Δx_{i}} , Δy_{i} ≐ y_{i+1} − y_{i} , Δx_{i} ≐ x_{i+1} − x_{i} .  (1.6)

Table 1.2: Linear interpolation in C

#include<a s s e r t . h>
double l i n t e r p ( int n , double∗ x , double∗ y , double z ){
a s s e r t ( z>=x [ 0 ] && z<=x [ n−1] && n>1);
int i =0, j=n−1; /∗ binary search : ∗/
while ( j−i >1){ int m=( i+j ) / 2 ; i f ( z>x [m] ) i=m; else j=m; }
return y [ i ]+(y [ i +1]−y [ i ] ) / ( x [ i +1]−x [ i ] ) ∗ ( z−x [ i ] ) ;

}

Note that the search of the interval [x_{i} ≤ x ≤ x_{i+1}] in an ordered array {x_{i}} should be done with the binary search algoritm (also called half-interval search): the point x is compared to the middle element of the array, if it is less than the middle element, the algorithm repeats its action on the subarray to the left of the middle element, if it is greater, on the sub-array to the right. When the remaining sub-array is reduced to two elements, the interval is found. The average number of operations for a binary search is O(log n).

1.3.2 Quadratic spline

Quadratic spline is made of second order polynomials, conveniently written in the form

S_{i}(x) = y_{i} + p_{i}(x − x_{i}) + c_{i}(x − x_{i})(x − x_{i+1}) |_{i=1,...,n−1} ,  (1.7)

which identically satisfies the continuity conditions

S_{i}(x_{i}) = y_{i} , S_{i}(x_{i+1}) = y_{i+1} |_{i=1,...,n−1} .  (1.8)

Substituting (1.7) into the derivative continuity condition,

S′_{i}(x_{i+1}) = S′_{i+1}(x_{i+1}) |_{i=1,...,n−2} ,  (1.9)

gives n − 2 equations for n − 1 unknown coefficients c_{i},

p_{i} + c_{i}Δx_{i} = p_{i+1} − c_{i+1}Δx_{i+1} |_{i=1,...,n−2} .  (1.10)

One coefficient can be chosen arbitrarily, for example c_{1} = 0. The other coefficients can now be calculated recursively,

c_{i+1} = \frac{1}{Δx_{i+1}} (p_{i+1} − p_{i} − c_{i}Δx_{i}) |_{i=1,...,n−2} .  (1.11)

Alternatively, one can choose c_{n−1} = 0 and make the backward-recursion

c_{i} = \frac{1}{Δx_{i}} (p_{i+1} − p_{i} − c_{i+1}Δx_{i+1}) |_{i=n−2,...,1} .  (1.12)

In practice, unless you know what your c_{1} (or c_{n−1}) is, it is better to run both recursions and then average the resulting c’s. This amounts to first running the forward-recursion from c_{1} = 0 and then the backward recursion from \frac{1}{2}c_{n−1}.

The optimized form (1.7) of the quadratic spline can aslo be written in the ordinary form, suitable for differentiation and integration

S_{i}(x) = y_{i} + b_{i}(x − x_{i}) + c_{i}(x − x_{i})^{2} , where b_{i} = p_{i} − c_{i}Δx_{i} . (1.13)

An implementation of quadratic spline in C is listed in Table 1.3.2

1.3.3 Cubic spline

Cubic splines are made of third order polynomials,

S_{i}(x) = y_{i} + b_{i}(x − x_{i}) + c_{i}(x − x_{i})^{2} + d_{i}(x − x_{i})^{3} .  (1.14)

1.3. SPLINE INTERPOLATION 5

Table 1.3: Quadratic spline in C

#include <s t d l i b . h>
#include <a s s e r t . h>
typedef struct { int n ; double ∗x , ∗y , ∗b , ∗ c ; } q s p l i n e ;
q s p l i n e ∗ q s p l i n e a l l o c ( int n , double∗ x , double∗ y ){
// b u i l d s q u a d r a t i c s p l i n e
q s p l i n e ∗ s = ( q s p l i n e ∗) malloc ( sizeof ( q s p l i n e ) ) ; // s p l i n e
s−>b = ( double ∗) malloc ( ( n−1)∗ sizeof ( double ) ) ; // b i
s−>c = ( double ∗) malloc ( ( n−1)∗ sizeof ( double ) ) ; // c i
s−>x = ( double ∗) malloc ( n∗ sizeof ( double ) ) ; // copy o f x i
s−>y = ( double ∗) malloc ( n∗ sizeof ( double ) ) ; // copy o f y i
s−>n = n ; for ( int i =0; i<n ; i ++){s−>x [ i ]=x [ i ] ; s−>y [ i ]=y [ i ] ; }
int i ; double p [ n−1] , h [ n −1]; //VLA from C99
for ( i =0; i<n−1; i ++){h [ i ]=x [ i +1]−x [ i ] ; p [ i ]=(y [ i +1]−y [ i ] ) / h [ i ] ; }
s−>c [ 0 ] = 0 ; // r e c u r s i o n up :
for ( i =0; i<n−2; i++)s−>c [ i +1]=(p [ i +1]−p [ i ]−s−>c [ i ] ∗ h [ i ] ) / h [ i +1];
s−>c [ n−2]/=2; // r e c u r s i o n down :
for ( i=n−3; i >=0;i −−)s−>c [ i ]=(p [ i +1]−p [ i ]−s−>c [ i +1]∗h [ i +1])/h [ i ] ;
for ( i =0; i<n−1; i++)s−>b [ i ]=p [ i ]−s−>c [ i ] ∗ h [ i ] ;
return s ; }
double q s p l i n e e v a l ( q s p l i n e ∗s , double z ){ // e v a l u a t e s s ( z )
a s s e r t ( z>=s−>x [ 0 ] && z<=s−>x [ s−>n − 1 ] ) ;
int i =0, j=s−>n−1; // binary search :
while ( j−i >1){ int m=( i+j ) / 2 ; i f ( z>s−>x [m] ) i=m; else j=m; }
double h=z−s−>x [ i ] ;
return s−>y [ i ]+h∗( s−>b [ i ]+h∗s−>c [ i ] ) ; }// i n e r p o l a t i n g polynomial
void q s p l i n e f r e e ( q s p l i n e ∗ s ){ // f r e e the a l l o c a t e d memory
f r e e ( s−>x ) ; f r e e ( s−>y ) ; f r e e ( s−>b ) ; f r e e ( s−>c ) ; f r e e ( s ) ; }

This form automatically satisfies the first half of continuity conditions (1.3): S_{i}(x_{i}) = y_{i}. The second half of continuity conditions (1.3), S_{i}(x_{i+1}) = y_{i+1}, and the continuity of the first and second derivatives (1.4) give a set of equations,

y_{i} + b_{i}h_{i} + c_{i}h_{i}^{2} + d_{i}h_{i}^{3} = y_{i+1} , i = 1, …, n − 1
b_{i} + 2c_{i}h_{i} + 3d_{i}h_{i}^{2} = b_{i+1} , i = 1, …, n − 2
2c_{i} + 6d_{i}h_{i} = 2c_{i+1} , i = 1, …, n − 2  (1.15)

where

h_{i} = x_{i+1} − x_{i} .  (1.16)

The set of equations (1.15) is a set of 3n − 5 linear equations for the 3(n − 1) unknown coefficients {a_{i}, b_{i}, c_{i} | i = 1, …, n − 1}. Therefore two more equations should be added to the set to find the coefficients. If the two extra equations are also linear, the total system is linear and can be easily solved.

The spline is called natural if the extra conditions are given as vanishing second derivative at the end-points,

S′′(x_{1}) = S′′(x_{n}) = 0 ,  (1.17)

which gives

c_{1} = 0 ,
c_{n−1} + 3d_{n−1}h_{n−1} = 0 .  (1.18)

Solving the first two equations in (1.15) for c_{i} and d_{i} gives^{1}

c_{i}h_{i} = −2b_{i} − b_{i+1} + 3p_{i} ,
d_{i}h_{i}^{2} = b_{i} + b_{i+1} − 2p_{i} ,  (1.19)

where p_{i} ≐ \frac{Δy_{i}}{h_{i}}. The natural conditions (1.18) and the third equation in (1.15) then produce the following tridiagonal system of n linear equations for the n coefficients b_{i},

2b_{1} + b_{2} = 3p_{1} ,
b_{i} + (2\frac{h_{i}}{h_{i+1}} + 2)b_{i+1} + \frac{h_{i}}{h_{i+1}} b_{i+2} = 3(p_{i} + p_{i+1} \frac{h_{i}}{h_{i+1}} ) |_{i=1,...,n−2} ,
b_{n−1} + 2b_{n} = 3p_{n−1} ,  (1.20)

or, in the matrix form,

( \begin{matrix}D_{1} & Q_{1} & 0 & 0 & … \\ 1 & D_{2} & Q_{2} & 0 & … \\ 0 & 1 & D_{3} & Q_{3} & … \\ ⋮ & ⋮ & ... & ... & ... \\ … & … & 0 & 1 & D_{n} \end{matrix} ) ( b_{1} ; ⋮ ; ⋮ ; b_{n} ) = ( B_{1} ; ⋮ ; ⋮ ; B_{n} )  (1.21)

where the elements D_{i} at the main diagonal are

D_{1} = 2 , D_{i+1} = 2\frac{h_{i}}{h_{i+1}} + 2 |_{i=1,...,n−2} , D_{n} = 2 ,  (1.22)

the elements Q_{i} at the above-main diagonal are

Q_{1} = 1 , Q_{i+1} = \frac{h_{i}}{h_{i+1}} |_{i=1,...,n−2} ,  (1.23)

and the right-hand side terms B_{i} are

B_{1} = 3p_{1} , B_{i+1} = 3 (p_{i} + p_{i+1} \frac{h_{i}}{h_{i+1}} ) |_{i=1,...,n−2} , B_{n} = 3p_{n−1} . (1.24)

1introducing an auxiliary coefficient b_{n}.

This system can be solved by one run of Gauss elimination and then a run of back-substitution. After a run of Gaussian elimination the system becomes

( \begin{matrix}D̃_{1} & Q_{1} & 0 & 0 & … \\ 0 & D̃_{2} & Q_{2} & 0 & … \\ 0 & 0 & D̃_{3} & Q_{3} & … \\ ⋮ & ⋮ & ... & ... & ... \\ … & … & 0 & 0 & D̃_{n} \end{matrix} ) ( b_{1} ; ⋮ ; ⋮ ; b_{n} ) = ( B̃_{1} ; ⋮ ; ⋮ ; B̃_{n} ) ,  (1.25)

where

D̃_{1} = D_{1} , D̃_{i} = D_{i} − Q_{i−1}/D̃_{i−1} |_{i=2,...,n} ,  (1.26)

and

B̃_{1} = B_{1} , B̃_{i} = B_{i} − B̃_{i−1}/D̃_{i−1} |_{i=2,...,n} .  (1.27)

The triangular system (1.25) can be solved by a run of back-substitution,

b_{n} = \frac{1}{D̃_{n}} B̃_{n} , b_{i} = \frac{1}{D̃_{i}} (B̃_{i} − Q_{i}b_{i+1}) |_{i=n−1,...,1} .  (1.28)

An implementation of quadratic spline in C is listed in Table 1.3.3

1.3.4 Akima sub-spline interpolation

Akima sub-spline is an interpolating function in the form of a piecewise cubic polynomial [1], similar to a cubic spline,

A(x)|_{x∈[x_{i},x_{i+1}]} = a_{i} + b_{i}(x − x_{i}) + c_{i}(x − x_{i})^{2} + d_{i}(x − x_{i})^{3} ≐ A_{i}(x) . (1.29)

However, unlike the cubic spline, Akima sub-spline dispenses with the demand of maximal differentiability of the spline—in this case, the continuity of the second derivative—hence the name sub-spline. Instead of achieving maximal differentiability Akima sub-splines try to reduce the wiggling the ordinary splines are typically prone to (see figure 1.2).

First let us note that the coefficients {a_{i}, b_{i}, c_{i}, d_{i}} in eq. (1.29) are determined by the values of the derivatives A′_{i} ≐ A′(x_{i}) of the sub-spline through the continuity conditions for the sub-spline and its first derivative,

A_{i}(x_{i}) = y_{i}, A′_{i}(x_{i}) = A′_{i}, A_{i}(x_{i+1}) = y_{i+1}, A′_{i}(x_{i+1}) = A′_{i+1}.  (1.30)

Inserting (1.29) into (1.30) and solving for the coefficients gives

a_{i} = y_{i}, b_{i} = A′_{i}, c_{i} = \frac{3p_{i} − 2A′_{i} − A′_{i+1}}{Δx_{i}} , d_{i} = \frac{A′_{i} + A′_{i+1} − 2p_{i}}{(Δx_{i})^{2}} , (1.31)

where p_{i} ≐ Δy_{i}/Δx_{i}, Δy_{i} ≐ y_{i+1} − y_{i}, Δx_{i} ≐ x_{i+1} − x_{i}.

Figure 1.2: A cubic spline (solid line) showing the typical wiggles, compared to the Akima sub-spline (dashed line) where the wiggles are essentially removed.

In the ordinary qubic spline the derivatives A′_{i} are determined by the continuity condition of the second derivative of the spline. Sub-splines do without this continuity condition and can instead use the derivatives as free parameters to be chosen to satisfy some other condition.

Akima suggested to minimize the wiggling by choosing the derivatives as linear combinations of the nearest slopes,

A′_{i} = \frac{w_{i+1}p_{i−1} + w_{i−1}p_{i}}{w_{i+1} + w_{i−1}} , if w_{i+1} + w_{i−1} ≠ 0 ,  (1.32)
A′_{i} = \frac{p_{i−1} + p_{i}}{2} , if w_{i+1} + w_{i−1} = 0 ,  (1.33)

where the weights w_{i} are given as

w_{i} = |p_{i} − p_{i−1}| .  (1.34)

The idea is that if three points lie close to a line, the sub-spline in this vicinity has to be close to this line. In other words, if |p_{i} − p_{i−1}| is small, the nearby derivatives must be close to p_{i}.

The first two and the last two points need a special prescription, for example (naively) one can simply use

A′_{1} = p_{1}, A′_{1} = \frac{1}{2}p_{1} + \frac{1}{2}p2, A′_{n} = p_{n−1}, A′_{n−1} = \frac{1}{2}p_{n−1} + \frac{1}{2}p_{n−2}. (1.35)

Table (1.5) shows a javascript implementation of this algorithm.

1.4 Other forms of interpolation

Other forms of interpolation can be constructed by choosing different classes of interpolating functions, for example, rational function interpolation, trigonometric interpolation, wavelet interpolation etc.

Sometimes not only the values of the function are tabulated but also the values of its derivative. This extra information can be taken advantage of when constructing the interpolation function.

1.5 Multivariate interpolation

Interpolation of a function in more than one varibale is called multivariate interpolation. The function of interest is represented as a set of discrete points in a multidimentional space. The points may or may not lie on a regular grid.

1.5.1 Nearest-neighbor interpolation

Nearest-neighbor interpolation approximates the value of the function at a non-tabulated point by the value at the nearest tabulated point, yielding a piecewise-constant interpolating function. It can be used for both regular and irregular grids.

1.5.2 Piecewise-linear interpolation

Piecewise-linear interpolation is used to interpolate functions of two variables tabulated on irregular grids. The tabulated 2D region is triangulated – subdivided into a set of non-intersecting triangles whose union is the original region. Inside each triangle the interpolating function S(x, y) is taken in the linear form,

S(x, y) = a + bx + cy ,  (1.36)

where the three constants are determined by the three conditions that the interpolating function is equal the tabulated values at the three vertices of the triangle.

1.5.3 Bi-linear interpolation

Bi-linear interpolation is used to interpolate functions of two variables tabulated on regular rectilinear 2D grids.

The interpolating function B(x, y) inside each of the grid rectangles is taken as a product of two linear functions of x and y respectively,

B(x, y) = (α + βx)(γ + δy) ≐ a + bx + cy + dxy ,  (1.37)

where the four constants a, b, c, d are obtained from the four conditions that the interpolating function is equal the tabulated values at the four nearest tabulated points (which are the vertices of the given grid rectangle).

1.5. MULTIVARIATE INTERPOLATION 11

Table 1.4: Cubic spline in C

#include<s t d l i b . h>
#include<a s s e r t . h>
typedef struct { int n ; double ∗x , ∗ y , ∗ b , ∗ c , ∗ d ; } c s p l i n e ;
c s p l i n e ∗ c s p l i n e a l l o c ( int n , double ∗x , double ∗y )
{// b u i l d s n a t u r a l c u b i c s p l i n e
c s p l i n e ∗ s = ( c s p l i n e ∗) malloc ( sizeof ( c s p l i n e ) ) ;
s−>x = ( double ∗) malloc ( n∗ sizeof ( double ) ) ;
s−>y = ( double ∗) malloc ( n∗ sizeof ( double ) ) ;
s−>b = ( double ∗) malloc ( n∗ sizeof ( double ) ) ;
s−>c = ( double ∗) malloc ( ( n−1)∗ sizeof ( double ) ) ;
s−>d = ( double ∗) malloc ( ( n−1)∗ sizeof ( double ) ) ;
s−>n = n ; for ( int i =0; i<n ; i ++){s−>x [ i ]=x [ i ] ; s−>y [ i ]=y [ i ] ; }
double h [ n−1] ,p [ n −1]; // VLA
for ( int i =0; i<n−1; i++) h [ i ]=x [ i +1]−x [ i ] ;
for ( int i =0; i<n−1; i++) p [ i ]=(y [ i +1]−y [ i ] ) / h [ i ] ;
double D[ n ] , Q[ n−1] , B[ n ] ; // b u i l d i n g the t r i d i a g o n a l system :
D[ 0 ] = 2 ; for ( int i =0; i<n−2; i++)D[ i +1]=2∗h [ i ] / h [ i +1]+2; D[ n−1]=2;
Q[ 0 ] = 1 ; for ( int i =0; i<n−2; i++)Q[ i +1]=h [ i ] / h [ i +1];
for ( int i =0; i<n−2; i++)B[ i +1]=3∗(p [ i ]+p [ i +1]∗h [ i ] / h [ i + 1 ] ) ;
B[0]=3∗ p [ 0 ] ; B[ n−1]=3∗p [ n −2]; //Gauss e l i m i n a t i o n :
for ( int i =1; i<n ; i ++){ D[ i ]−=Q[ i −1]/D[ i −1]; B[ i ]−=B[ i −1]/D[ i −1]; }
s−>b [ n−1]=B[ n−1]/D[ n −1]; // back−s u b s t i t u t i o n :
for ( int i=n−2; i >=0;i −−) s−>b [ i ]=(B[ i ]−Q[ i ] ∗ s−>b [ i +1])/D[ i ] ;
for ( int i =0; i<n−1; i ++){
s−>c [ i ]=(−2∗s−>b [ i ]−s−>b [ i +1]+3∗p [ i ] ) / h [ i ] ;
s−>d [ i ]=( s−>b [ i ]+s−>b [ i +1]−2∗p [ i ] ) / h [ i ] / h [ i ] ;

}

return s ;

}

double c s p l i n e e v a l ( c s p l i n e ∗s , double z ){ // e v a l u a t i o n o f s p l i n e
a s s e r t ( z>=s−>x [ 0 ] && z<=s−>x [ s−>n − 1 ] ) ;
int i =0, j=s−>n−1;// binary search f o r the i n t e r v a l f o r z :
while ( j−i >1){ int m=( i+j ) / 2 ; i f ( z>s−>x [m] ) i=m; else j=m; }
double h=z−s−>x [ i ] ; // c a l c u l a t e the i n e r p o l a t i n g s p l i n e :
return s−>y [ i ]+h∗( s−>b [ i ]+h∗( s−>c [ i ]+h∗s−>d [ i ] ) ) ;

}

void c s p l i n e f r e e ( c s p l i n e ∗ s ){ // f r e e the a l l o c a t e d memory
f r e e ( s−>x ) ; f r e e ( s−>y ) ; f r e e ( s−>b ) ; f r e e ( s−>c ) ; f r e e ( s−>d ) ; f r e e ( s ) ; }

12 CHAPTER 1. INTERPOLATION

Table 1.5: Javascript implementation of Akima sub-spline.

Array . prototype . i t e r a t o r =f u n c t i o n ( )
{ for ( var i =0; i<t h i s . length ; i++)y i e l d i }
f u n c t i o n akima ( x , y ){
var n=x . length ;
var h=[(x [ i +1]−x [ i ] ) for ( i in x ) i f ( i<n − 1 ) ] ;
var p=[(y [ i +1]−y [ i ] ) / h [ i ] for ( i in h ) ] ;
var A=new Array ( n ) , c=new Array (n−1) , d=new Array (n−1);
A[0]= p [ 0 ] ; A[ 1 ] = ( p [0]+ p [ 1 ] ) / 2 ;
A[ n−1]=p [ n −2]; A[ n−2]=(p [ n−2]+p [ n −3])/2;
for ( var i =2; i<n−2; i++)
{ w1=Math . abs ( p [ i +1]−p [ i ] ) ; w2=Math . abs ( p [ i −1]−p [ i − 2 ] ) ;
i f (w1+w2==0) A[ i ]=(p [ i −1]+p [ i ] ) / 2 ;
else A[ i ]=(w1∗p [ i −1]+w2∗p [ i ] ) / ( w1+w2 ) ;

}

for ( var i =0; i<n−1; i++)
{ c [ i ]=(3∗p [ i ]−2∗A[ i ]−A[ i +1])/h [ i ] ;
d [ i ]=(A[ i +1]+A[ i ]−2∗p [ i ] ) / h [ i ] / h [ i ] ;

}

var s p l i n e=f u n c t i o n ( z )
{ i f ( z<x [ 0 ] | | z>x [ n−1]) throw ” s p l i n e : out o f range ” ;
var i =0, j=n−1;
while ( j−i >1){
var mid=Math . round ( ( i+j )/2 ) ;
i f ( z>x [ mid ] ) i=mid ; else j=mid ; }
var dx=z−x [ i ] ;
return y [ i ]+dx ∗(A[ i ]+dx ∗( c [ i ]+dx∗d [ i ] ) ) ;

}

return s p l i n e ;

}

Numerical integration

2.1 Introduction

Numerical integration constitutes a broad family of algorithms to compute a numerical approximation to a definite (Riemann) integral.

Generally, the integral is approximated by a weighted sum of function values within the domain of integration,

∫_{a}^{b} f(x)dx ≈ ∑_{i=1}^{n} w_{i}f(x_{i}) .  (2.1)

Expression (2.1) is often referred to as quadrature (cubature for multidimensional integrals) or rule. The abscissas x_{i} (also called nodes) and the weights w_{i} of a quadrature are usually optimized—using one of a large number of different strategies—to suit a particular class of integration problems.

The best quadrature algorithm for a given problem depends on several factors, in particular on the integrand. Different classes of integrands generally require different quadratures for the most effective calculation.

A popular numerical integration library is QUADPACK [14]. It includes general purpose routines—like QAGS, basen on an adaptive Gauss–Kronrod quadrature with acceleration—as well as a number of specialized routines. The GNU scientific library [4] (GSL) implements most of the QUADPACK routines and in addition includes a modern general-purpose adaptive routine CQUAD based on Clenshaw-Curtis quadratures [9].

In the following we shall consider some of the popular numerical integration algorithms.

2.2 Rectangle and trapezium rules

In mathematics, the Reimann integral is generally defined in terms of Riemann sums [15]. If the integration interval [a, b] is partitioned into n subintervals,

a = t_{0} < t_{1} < t_{2} < ⋯ < t_{n} = b .  (2.2)

the Riemann sum is defined as

∑_{i=1}^{n} f(x_{i})Δx_{i} ,  (2.3)

where x_{i} ∈ [t_{i−1}, t_{i}] and Δx_{i} = t_{i} − t_{i−1}. Geometrically a Riemann sum can be interpreted as the area of a collection of adjucent rectangles with widths Δx_{i} and heights f(x_{i}).

The Rieman integral is defined as the limit of a Riemann sum as the mesh—the length of the largest subinterval—of the partition approaches zero. Specifically, the number denoted as

∫_{a}^{b} f(x)dx  (2.4)

is called the Riemann integral, if for any ϵ > 0 there exists δ > 0 such that for any partition (2.2) with max Δx_{i} < δ we have

|∑_{i=1}^{n} f(x_{i})Δx_{i} − ∫_{a}^{b} f(x)dx| < ϵ .  (2.5)

A definite integral can be interpreted as the net signed area bounded by the graph of the integrand.

Now, the n-point rectangle quadrature is simply the Riemann sum (2.3),

∫_{a}^{b} f(x)dx ≈ ∑_{i=1}^{n} f(x_{i})Δx_{i} ,  (2.6)

where the node x_{i} is often (but not always) taken at the middle of the corresponding subinterval, x_{i} = t_{i−1} + \frac{1}{2}Δx_{i}, and the subintervals are often (but not always) chosen equal, Δx_{i} = (b − a)/n. Geometrically the n-point rectangle rule is an approximation to the integral given by the area of a collection of n adjucent equal rectangles whose heights are determined by the values of the function (at the middle of the rectangle).

An n-point trapezium rule uses instead a collection of trapezia fitted under the graph,

∫_{a}^{b} f(x)dx ≈ ∑_{i=1}^{n} \frac{f(t_{i−1}) + f(t_{i})}{2}Δx_{i} .  (2.7)

Importantly, the trapezium rule is the average of two Riemann sums,

∑_{i=1}^{n} \frac{f(t_{i−1}) + f(t_{i})}{2}Δx_{i} = \frac{1}{2} ∑_{i=1}^{n} f(t_{i−1})Δx_{i} + \frac{1}{2} ∑_{i=1}^{n} f(t_{i})Δx_{i} .  (2.8)

2.3. QUADRATURES WITH REGULARLY SPACED ABSCISSAS 15

Rectangle and trapezium quadratures both have the important feature of closely following the very mathematical definition of the integral as the limit of the Riemann sums. Therefore—disregarding the round-off errors—these two rules cannot fail if the integral exists.

For certain partitions of the interval the rectangle and trapezium rules coincide. For example, for the nodes

x_{i} = a + (b − a)\frac{i − \frac{1}{2}}{n} , i = 1, …, n  (2.9)

both rules give the same quadrature with equal weights, w_{i} = (b − a)/n,

∫_{a}^{b} f(x)dx ≈ \frac{b − a}{n} ∑_{i=1}^{n} f (a + (b − a)\frac{i − \frac{1}{2}}{n} ) .  (2.10)

Rectangle and trapezium quadratures are rarely used on their own— because of the slow convergence—but they often serve as the basis for more advanced quadratures, for example adaptive quadratures and variable transformation quadratures considered below.

2.3 Quadratures with regularly spaced abscissas

A quadrature (2.1) with n predefined nodes x_{i} has n free parameters: the weights w_{i}. A set of n parameters can generally be tuned to satisfy n conditions. The archetypal set of conditions in quadratures is that the quadrature integrates exactly a set of n functions,

{ϕ_{1}(x), …, ϕ_{n}(x)} .  (2.11)

This leads to a set of n equations,

∑_{i=1}^{n} w_{i}ϕ_{k}(x_{i}) = I_{k} |_{k=1,...,n} ,  (2.12)

where the integrals

I_{k} ≐ ∫_{a}^{b} ϕ_{k}(x)dx  (2.13)

are assumed to be known. Equations (2.12) are linear in w_{i} and can be easily solved.

Since integration is a linear operation, the quadrature will then also integrate exactly any linear combination of functions (2.11).

A popular choice for predefined nodes is a closed set—that is, including the end-points of the interval—of evenly spaced abscissas,

x_{i} = a + \frac{i − 1}{n − 1}(b − a) |_{i=1,...,n} .  (2.14)

However, in practice it often happens that the integrand has an integrable sinfularity at one or both ends of the interval. In this case one can choose an open set of equidistant nodes,

x_{i} = a + \frac{i − \frac{1}{2}}{n} (b − a) |_{i=1,...,n} .  (2.15)

The set of functions to be integrated exactly is generally chosen to suite the properties of the integrands at hand: the integrands must be well represented by linear combinations of the chosen functions.

2.3.1 Classical quadratures

Suppose the integrand can be well represented by the first few terms of its Taylor series,

f(x) = ∑_{k=0}^{∞} \frac{f^{(k)}(a)}{k!}(x − a)^{k} ,  (2.16)

where f^{(k)} is the k-th derivative of the integrand. This is often the case for analytic—that is, infinitely differentiable—functions. For such integrands one can obviously choose polynomials

{1, x, x^{2}, …, x^{n−1}}  (2.17)

as the set of functions to be integrated exactly.

This leads to the so called classical quadratures: quadratures with regularly spaced abscissas and polynomials as exactly integrable functions.

An n-point classical quadrature integrates exactly the first n terms of the function’s Taylor expansion (2.16). The x^{n} order term will not be integrated exactly and will lead to an error of the quadrature. Thus the error E_{n} of the n-point classical quadrature is on the order of the integral of the x^{n} term in (2.16),

E_{n} ≈ ∫_{a}^{b} \frac{f^{(n)}(a)}{n!}(x − a)^{n}dx = \frac{f^{(n)}(a)}{(n + 1)!}h^{n+1} ∝ h^{n+1} ,  (2.18)

where h = b − a is the length of the integration interval. A quadrature with the error of the order h^{n+1} is often called a degree-n quadrature.

If the integrand is smooth enough and the length h is small enough a classical quadrature with not so large n can provide a good approximation for the integral. However, for large n the weights of classical quadratures tend to have alternating signs, which leads to large round-off errors, which in turn negates the potentially higher accuracy of the quadrature. Again, if the integrand violates the assumption of Taylor expansion—for example by having an integrable singularity inside the integration interval—the higher order quadratures may perform poorly.

Table 2.1: Maxima script to calculate analytically the weights of an n-point classical quadrature with predefined abscissas in the interval [0, 1].

n: 8; xs: makelist((i-1)/(n-1),i,1,n); /* nodes: adapt to your needs */
ws: makelist(concat(w,i),i,1,n);
ps: append([1],makelist(x^i,i,1,n-1)); /* polynomials */
fs: makelist(buildq([i:i,ps:ps],lambda([x],ps[i])),i,1,n);
integ01: lambda([f],integrate(f(x),x,0,1));
Is: maplist(integ01,fs); /* calculate the integrals */
eq: lambda([f],lreduce("+",maplist(f,xs)*ws));
eqs: maplist(eq,fs)-Is; /* build equations */
solve(eqs,ws); /* solve for the weights */

Classical quadratures are mostly of historical interest nowadays. Alternative methods—such as quadratures with optimized abscissas, adaptive, and variable transformation quadratures—are more stable and accurate and are normally preferred to classical quadratures.

Classical quadratures with equally spaced abscissas—both closed and open sets—are generally referred to as Newton-Cotes quadratures. An interested reader can generate Newton—Cotes quadratures of any degree n using the Maxima script in Table (2.1).

2.4 Quadratures with optimized abscissas

In quadratures with optimized abscissas not only the weights w_{i} but also the abscissas x_{i} are chosen optimally. The number of free parameters is thus 2n and one can choose a set of 2n functions,

{ϕ_{1}(x), …, ϕ_{2n}(x)} ,  (2.19)

to be integrated exactly. This gives a system of 2n equations, linear in w_{i} and non-linear in x_{i},

∑_{i=1}^{n} w_{i}ϕ_{k}(x_{i}) = I_{k} |_{k=1,...,2n} ,  (2.20)

where again

I_{k} ≐ ∫_{a}^{b} ϕ_{k}(x)dx .  (2.21)

The weights and abscissas of the quadrature can be determined by solving this system of equations^{1}.

1Here is, for example, an n = 2 quadrature with optimized abscissas,

∫_{−1}^{1} f(x)dx ≈ f (−\sqrt{\frac{1}{3}} ) + f (+\sqrt{\frac{1}{3}} ) .  (2.22)

Although quadratures with optimized abcissas are generally of much higher order, 2n−1 compared to n−1 for non-optimal abscissas, the optimal points generally can not be reused at the next iteration in an adaptive algorithm.

2.4.1 Gauss quadratures

Gauss quadratures deal with a slightly more general form of integrals,

∫_{a}^{b} ω(x)f(x)dx ,  (2.23)

where ω(x) is a positive weight function. For ω(x) = 1 the problem is the same as considered above. Popular choices of the weight function include ω(x) = (1−x^{2})^{±1/2}, exp(−x), exp(−x^{2}) and others. The idea is to represent the integrand as a product ω(x)f(x) such that all the difficulties go into the weight function ω(x) while the remaining factor f(x) is smooth and well represented by polynomials.

An N-point Gauss quadrature is a quadrature with optimized abcissas,

∫_{a}^{b} ω(x)f(x)dx ≈ ∑_{i=1}^{N} w_{i}f(x_{i}) ,  (2.24)

which integrates exactly a set of 2N polynomials of the orders 1, …, 2N − 1 with the given weight ω(x).

Fundamental theorem

There is a theorem stating that there exists a set of polynomials p_{n}(x), orthogonal on the interval [a, b] with the weight function ω(x),

∫_{a}^{b} ω(x)p_{n}(x)p_{k}(x) ∝ δ_{nk} .  (2.25)

Now, one can prove that the optimal nodes for the N-point Gauss quadrature are the roots of the polynomial p_{N}(x),

p_{N}(x_{i}) = 0 .  (2.26)

The idea behind the proof is to consider the integral

∫_{a}^{b} ω(x)q(x)p_{N}(x)dx = 0 ,  (2.27)

where q(x) is an arbitrary polynomial of degree less than N. The quadrature should represent this integral exactly,

∑_{i=1}^{N} q(x_{i})p_{N}(x_{i}) = 0 .  (2.28)

Apparently this is only possible if x_{i} are the roots of p_{N} ■.

Calculation of nodes and weights

A neat algorithm—usually refered to as Golub-Welsch [8] algorithm—for calculation of the nodes and weights of a Gauss quadrature is based on the symmetric form of the three-term reccurence relation for orthogonal polynomials,

xp_{n−1}(x) = β_{n}p_{n}(x) + α_{n}p_{n−1}(x) + β_{n−1}p_{n−2}(x) ,  (2.29)

where p_{−1}(x) ≐ 0, p_{1}(x) ≐ 1, and n = 1, …, N. This reccurence relation can be written in the matrix form,

xp(x) = Jp(x) + β_{N}p_{N}(x)e_{N} ,  (2.30)

where p(x) ≐ {p_{0}(x), …, p_{N−1}(x)}^{T}, e_{N} = {0, …, 0, 1}^{T}, and the tridiagonal matrix J — usually refered to as Jacobi matrix or Jacobi operator — is given as

J = ( \begin{matrix}α_{1} & β_{1} \\ β_{1} & α_{2} & β_{2} \\ & β_{2} & α_{3} & β_{3} \\ & & ... & ... \\ & & & β_{N−1} & α_{N} \end{matrix} ) .  (2.31)

Substituting the roots x_{i} of p_{N}—that is, the set {x_{i} | p_{N}(x_{i}) = 0}— into the matrix equation (2.30) leads to eigenvalue problem for the Jacobi matrix,

Jp(x_{i}) = x_{i}p(x_{i}) .  (2.32)

Thus, the nodes of an N-point Gauss quadrature—the roots of the polynomial p_{N}—are the eigenvalues of the Jacobi matrix J and can be calculated by a standard diagonalization^{2} routine ■.

The weights can be obtained considering N integrals,

∫_{a}^{b} ω(x)p_{n}(x)dx = δ_{n0} ∫_{a}^{b} ω(x)dx , n = 0, …, N − 1 .  (2.33)

2A symmetric tridiagonal matrix can be diagonalized very effectively using the QR/RL algorithm.

Applying our quadrature gives the matrix equation,

Pw = e_{1} ∫_{a}^{b} ω(x)dx ,  (2.34)

where w ≐ {w_{1}, …, w_{N}}^{T}, e_{1} = {1, 0, …, 0}^{T}, and

P ≐ ( \begin{matrix}p_{0}(x_{1}) & … & p_{0}(x_{N}) \\ p_{1}(x_{1}) & … & p_{1}(x_{N}) \\ … & … & … \\ p_{N−1}(x_{1}) & … & p_{N−1}(x_{N}) \end{matrix} ) .  (2.35)

Equation (2.34) is linear in w_{i} and can be solved directly. However, if diagonalization of the Jacobi matrix provided the normalized eigenvectors, the weigths can be readily obtained using the following method.

The matrix P apparently consists of non-normalized column eigenvectors of the matrix J. The eigenvectors are orthogonal and therefore P^{T} P is a diagonal matrix with positive elements. Multiplying (2.34) by P^{T} and then by (P^{T} P)^{−1} from the left gives

w = (P^{T} P)^{−1}P^{T} e_{1} ∫_{a}^{b} ω(x)dx .  (2.36)

From p_{0}(x) = 1 it follows that P^{T} e_{1} = {1, …, 1}^{T} and therefore

w_{i} = \frac{1}{(P^{T} P)_{ii}} ∫_{a}^{b} ω(x)dx .  (2.37)

Let the matrix V be the set of the normalized column eigenvectors of the matrix J. The matrix V is then connected with the matrix P through the normalization equation,

V = \sqrt{(P^{T} P)^{−1}}P .  (2.38)

Therefore, again taking into account that p_{0}(x) = 1, equation (2.37) can be written as

w_{i} = (V_{1i})^{2} ∫_{a}^{b} ω(x)dx ■.  (2.39)

Example: Gauss-Legendre quadrature

Gauss-Legendre quadrature deals with the weight ω(x) = 1 on the interval [−1, 1]. The associated polynomials are Legendre polynomials P_{n}(x), hence the name. Their reccurence relation is usually given as

(2n − 1)xP_{n−1}(x) = nP_{n}(x) + (n − 1)P_{n−2}(x) .  (2.40)

Table 2.2: An Octave function which calculates the nodes and weights of the N-point Gauss-Legendre quadrature and then integrates a given function.

function Q = g a u s s l e g e n d r e ( f , a , b ,N)
beta = . 5 . / sqrt (1 −(2∗(1:N−1)).ˆ ( −2)); % r e c c u r e n c e r e l a t i o n
J = diag ( beta , 1 ) + diag ( beta , −1);  % Jacobi matrix
[V,D] = eig ( J ) ;  % d i a g o n a l i z a t i o n o f J
x = diag (D) ; [ x , i ] = sort ( x ) ;  % s o r t e d nodes
w = V( 1 , i ) . ˆ 2 ∗ 2 ;  % weights
Q = w∗ f ( ( a+b)/2+(b−a )/2∗ x ) ∗ ( b−a ) / 2 ; % i n t e g r a l
endfunction ;

Rescaling the polynomials (preserving p_{0}(x) = 1) as

\sqrt{2n + 1}P_{n}(x) = p_{n}(x)  (2.41)

reduces this reccurence relation to the symmetric form (2.29),

xp_{n−1}(x) = \frac{1}{2} \frac{1}{\sqrt{1 − (2n)^{−2}}} p_{n}(x) + \frac{1}{2} \frac{1}{\sqrt{1 − (2(n − 1))^{−2}}} p_{n−2}(x) . (2.42)

Correspondingly, the coefficients in the matrix J are

α_{n} = 0 , \{β_{n} = \frac{1}{2} \frac{1}{\sqrt{1 − (2n)^{−2}}} | n = 1, …, N − 1\} .  (2.43)

The problem of finding the nodes and the weights of the N-point Gauss-Legendre quadrature is thus reduced to the eigenvalue problem for the Jacobi matrix with coefficients (2.43).

As an illustration of this algorithm Table (2.2) shows an Octave function which calculates the nodes and the weights of the N-point Gauss-Legendre quadrature and then integrates a given function.

2.4.2 Gauss-Kronrod quadratures

Generally, the error of a numerical integration is estimated by comparing the results from two rules of different orders. However, for ordinary Gauss quadratures the nodes for two rules of different orders almost never coinside. This means that one can not reuse the points of the lower order rule when calculating the hihger order rule.

Gauss-Kronrod algorithm [11] remedies this inefficiency. The points inherited from the lower order rule are reused in the higher order rule as predefined nodes (with n weights as free parameters), and then m more optimal points are added (m abscissas and m weights as free parameters). The order of the method is n + 2m − 1. The lower order rule becomes embedded—that is, it uses a subset of the nodes—into the higher order rule. On the next iteration the procedure is repeated.

Patterson [13] has tabulated nodes and weigths for several sequences of embedded Gauss-Kronrod rules.

2.5 Adaptive quadratures

Higher order quadratures suffer from round-off errors as the weights w_{i} generally have alternating signs. Again, using high order polynomials is dangerous as they typically oscillate wildly and may lead to Runge’s phenomenon. Therefore, if the error of the quadrature is yet too large for a quadrature with sufficiently large n, the best strategy is to subdivide the interval in two and then use the quadrature on the half-intervals. Indeed, if the error is of the order h^{k}, the subdivision would lead to reduced error, 2 (h/2)^{k} < h^{k}, if k > 1.

An adaptive quadrature is an algorithm where the integration interval is subdivided into adaptively refined subintervals until the given accuracy goal is reached.

Adaptive algorithms are usually built on pairs of quadrature rules – a higher order rule,

Q = ∑_{i} w_{i}f(x_{i}),  (2.44)

where w_{i} are the weights of the higher order rule and Q is the higher order estimate of the integral, and a lower order rule,

q = ∑_{i} v_{i}f(x_{i}),  (2.45)

where v_{i} are the weights of the lower order rule and q is the the lower order estimate of the integral. The difference between the higher order rule and the lower order rule gives an estimate of the error,

δQ = |Q − q| .  (2.46)

The integration result is accepted, if the error δQ is smaller than tolerance,

δQ < δ + ϵ|Q| ,  (2.47)

where δ is the absolute accuracy goal and ϵ is the relative accuracy goal of the integration.

If the error estimate is larger than tolerance, the interval is subdivided into two half-intervals and the procedure applies recursively to subintervals with the same relative accuracy goal ϵ and rescaled absolute accuracy goal δ/\sqrt{2}.

The points x_{i} are usually chosen such that the two quadratures use the same points, and that the points can be reused in the subsequent recursive steps. The reuse of the function evaluations made at the previous step of adaptive integration is very important for the efficiency of the algorithm. The equally-spaced abscissas naturally provide for such a reuse.

As an example, Table 2.3 shows an implementation of the described algorithm using

x_{i} = \{\frac{1}{6}, \frac{2}{6}, \frac{4}{6}, \frac{5}{6}\} (easily reusable points) ,  (2.48)
w_{i} = \{\frac{2}{6}, \frac{1}{6}, \frac{1}{6}, \frac{2}{6}\} (trapezium rule) ,  (2.49)
v_{i} = \{\frac{1}{4}, \frac{1}{4}, \frac{1}{4}, \frac{1}{4}\} (rectangle rule) .  (2.50)

During recursion the function values at the points #2 and #3 are inherited from the previous step and need not be recalculated.

The points and weights are cited for the rescaled integration interval [0, 1]. The transformation of the points and weights to the original interval [a, b] is given as

x_{i} → a + (b − a)x_{i} ,
w_{i} → (b − a)w_{i} .  (2.51)

This implementation calculates directly the Riemann sums and can therefore deal with integrable singularities, although rather inefficiently.

More efficient adaptive routines keep track of the subdivisions of the interval and the local errors [9]. This allows detection of singularities and switching in their vicinity to specifically tuned quadratures. It also allows better estimates of local and global errors.

2.6 Variable transformation quadratures

The idea behind variable transformation quadratures is to apply the given quadrature—either with optimimized or regularly spaced nodes—not to the original integral, but to a variable transformed integral [12],

∫_{a}^{b} f(x)dx = ∫_{t_{a}}^{t_{b}} f(g(t))g′(t)dt ≈ ∑_{i=1}^{N} w_{i}f(g(t_{i}))g′(t_{i}) ,  (2.52)

where the transformation x = g(t) is chosen such that the transformed integral better suits the given quadrature. Here g′ denotes the derivative and [t_{a}, t_{b}] is the corresponding interval in the new variable.

For example, the Gauss-Legendre quadrature assumes the integrand can be well represented with polynomials and performs poorly on integrals with integrable singularities like

I = ∫_{0}^{1} \frac{1}{2\sqrt{x}}dx .  (2.53)

24 CHAPTER 2. NUMERICAL INTEGRATION

Table 2.3: Recursive adaptive integrator in C

#include<math . h>
#include<a s s e r t . h>
#include<s t d i o . h>
double adapt24 ( double f ( double ) , double a , double b ,
double acc , double eps , double f2 , double f3 , int nrec )
{ a s s e r t ( nrec <1000000);
double f 1=f ( a+(b−a ) / 6 ) , f 4=f ( a+5∗(b−a ) / 6 ) ;
double Q=(2∗ f 1+f 2+f 3 +2∗ f 4 )/6∗( b−a ) , q=( f 1+f 4+f 2+f 3 )/4∗( b−a ) ;
double t o l e r a n c e=acc+eps ∗ f a b s (Q) , e r r o r=f a b s (Q−q ) ;
i f ( e r r o r < t o l e r a n c e ) return Q;
else {
double Q1=adapt24 ( f , a , ( a+b )/2 , acc / s q r t ( 2 . ) , eps , f1 , f2 , nrec +1);
double Q2=adapt24 ( f , ( a+b )/2 , b , acc / s q r t ( 2 . ) , eps , f3 , f4 , nrec +1);
return Q1+Q2 ; }

}

double adapt ( double f ( double ) , double a , double b ,
double acc , double eps )
{ double f 2=f ( a+2∗(b−a ) / 6 ) , f 3=f ( a+4∗(b−a ) / 6 ) ; int nrec =0;
return adapt24 ( f , a , b , acc , eps , f2 , f3 , nrec ) ;

}

int main ( ) // uses gcc nested f u n c t i o n s
{ int c a l l s =0; double a=0,b=1, acc =0.001 , eps =0.001;
double f ( double x ){ c a l l s ++; return 1/ s q r t ( x ) ; }; // nested f u n c t i o n
double Q=adapt ( f , a , b , acc , eps ) ; p r i n t f ( ”Q=%g c a l l s=%d\n” ,Q, c a l l s ) ;
return 0 ;

}

However, a simple varibale transformation x = t^{2} removes the singularity,

I = ∫_{0}^{1} dt ,  (2.54)

and the Gauss-Legendre quadrature for the transformed integral gives exact result. The price is that the transformed quadrature performs less effectively on smooth functions.

Some of the popular variable transformation quadratures are Clenshaw-Curtis [3], based on the transformation

∫_{−1}^{1} f(x)dx = ∫_{0}^{π} f(cos θ) sin θdθ ,  (2.55)

and “tanh-sinh” quadrature [12], based on the transformation

∫_{−1}^{1} f(x)dx = ∫_{−∞}^{∞} f (tanh (\frac{π}{2} sinh(t))) \frac{π}{2} \frac{cosh(t)}{cosh^{2} (\frac{π}{2} sinh(t))}dt . (2.56)

Generally, the equally spaced trapezium rule is used after the transformation.

2.7 Infinite intervals

One way to calculate an integral over infinite interval is to transform it by a variable sustitution into an integral over a finite interval. The latter can then be evaluated by ordinary integration methods. Table 2.4 lists several of such transformation.

Table 2.4: Variable transformations reducing infinite interval integrals into integrals over finite intervals.

∫_{−∞}^{+∞} f(x)dx = ∫_{−1}^{+1} f (\frac{t}{1 − t^{2}} ) \frac{1 + t^{2}}{(1 − t^{2})^{2}} dt ,  (2.57)
∫_{−∞}^{+∞} f(x)dx = ∫_{0}^{1} (f (\frac{1 − t}{t}) + f (−\frac{1 − t}{t})) \frac{dt}{t^{2}} , (2.58)
∫_{a}^{+∞} f(x)dx = ∫_{0}^{1} f (a + \frac{t}{1 − t}) \frac{1}{(1 − t)^{2}} dt ,  (2.59)
∫_{a}^{+∞} f(x)dx = ∫_{0}^{1} f (a + \frac{1 − t}{t}) \frac{dt}{t^{2}} ,  (2.60)
∫_{−∞}^{a} f(x)dx = ∫_{−1}^{0} f (a − \frac{t}{1 + t}) \frac{−1}{(1 + t)^{2}} dt ,  (2.61)
∫_{−∞}^{a} f(x)dx = ∫_{0}^{1} f (a − \frac{1 − t}{t}) \frac{dt}{t^{2}} .  (2.62)

26 CHAPTER 2. NUMERICAL INTEGRATION

Monte Carlo integration

3.1 Introduction

Monte Carlo integration is a quadrature (cubature) where the nodes are chosen randomly. Typically no assumptions are made about the smoothness of the integrand, not even that it is continuous.

Monte Carlo algorithms are particularly suited for multi-dimensional integrations where one of the problems is that the integration region, Ω, often has quite complicated boundary which can not be easily described by simple functions. However, it is usually much easier to find out whether a given point lies within the integration region or not. Therefore a popular strategy is to create an auxiliary rectangular volume, V, which encompasses the integration volume Ω, and an auxiliary function which coincides with the integrand inside the volume Ω and is equal zero outside. Then the integral of the auxiliary function over the auxiliary volume is equal the original integral.

However, the auxiliary function is generally non-continuous at the boundary; thus ordinary quadratures which assume continuity of the integrand are bound to have difficulties here. One the contrary the Monte-Carlo quadratures will do just as good (or as bad) as with continuous integrands.

A typical implementation of a Monte Carlo algorithm integrates the given function over a rectangular volume, specified by the coordinates of its ”lower-left” and ”upper-right” vertices, assuming the user has provided the encompassing volume with the auxiliary function.

Plain Monte Carlo algorithm distributes points uniformly throughout the integration region using uncorrelated pseudo-random sequences of points.

Adaptive algorithms, such as VEGAS and MISER, distribute points nonuniformly in an attempt to reduce integration error using correspondingly importance and stratified sampling.

Yet another strategy to reduce the error is to use correlated quasi-random sequences.

The GNU Scientific Library, GSL, implements a plain Monte Carlo integration algorithm; a stratified sampling algorithm, MISER; an importance sampling algorithm, VEGAS; and a number of quasi-random generators.

3.2 Plain Monte Carlo sampling

Plain Monte Carlo is a quadrature with random abscissas and equal weights,

∫_{V} f(x)dV ≈ w ∑_{i=1}^{N} f(x_{i}) ,  (3.1)

where x is a point in the multi-dimensional integration space. One free parameter, w, allows one condition to be satisfied: the quadrature must integrate exactly a constant function. This gives w = V/N,

∫_{V} f(x)dV ≈ \frac{V}{N} ∑_{i=1}^{N} f(x_{i}) ≐ V ⟨f⟩ .  (3.2)

Under the assumptions of the central limit theorem the error of the integration can be estimated as

error = V \frac{σ}{\sqrt{N}} ,  (3.3)

where σ is the variance of the sample,

σ^{2} = ⟨f^{2}⟩ − ⟨f⟩^{2} .  (3.4)

The familiar 1/\sqrt{N} convergence of a random walk process is quite slow: to reduce the error by a factor 10 requires 100-fold increase in the number of sample points.

Expression (3.3) provides only a statistical estimate of the error, which is not a strict bound; random sampling may not uncover all the important features of the function, resulting in an underestimate of the error.

A simple implementation of the plain Monte Carlo algorithm is shown in Table 3.1.

3.3 Importance sampling

Suppose the points are distributed not uniformly but with some density ρ(x) . That is, the number of points Δn in the volume ΔV around point x is given as

Δn = \frac{N}{V} ρ(x)ΔV,  (3.5)

where ρ is normalised such that ∫_{V} ρdV = V.

Table 3.1: Plain Monte Carlo integrator

#include <math . h>
#include <s t d l i b . h>
#define RND ( ( double ) rand ( ) /RAND MAX)
void randomx ( int dim , double ∗a , double ∗b , double ∗x )
{ for ( int i =0; i<dim ; i++) x [ i ]=a [ i ]+RND∗(b [ i ]−a [ i ] ) ; }
void plainmc ( int dim , double ∗a , double ∗b ,
double f ( double∗ x ) , int N, double∗ r e s u l t , double∗ e r r o r )
{ double V=1; for ( int i =0; i<dim ; i++) V∗=b [ i ]−a [ i ] ;
double sum=0, sum2=0, fx , x [ dim ] ;
for ( int i =0; i<N; i ++){ randomx ( dim , a , b , x ) ; fx=f ( x ) ;
sum+=fx ; sum2+=fx ∗ fx ; }
double avr = sum/N, var = sum2/N−avr ∗ avr ;
∗ r e s u l t = avr ∗V; ∗ e r r o r = s q r t ( var /N)∗V;

}

The estimate of the integral is then given as

∫_{V} f(x)dV ≈ ∑_{i=1}^{N} f(x_{i})ΔV_{i} = ∑_{i=1}^{N} f(x_{i})\frac{V}{Nρ(x_{i})} = V ⟨\frac{f}{ρ} ⟩ ,  (3.6)

where

ΔV_{i} = \frac{V}{Nρ(x_{i})}  (3.7)

is the volume-per-point at the point x_{i}.

The corresponding variance is now given by

σ^{2} = ⟨(\frac{f}{ρ} )^{2}⟩ − ⟨\frac{f}{ρ} ⟩^{2} .  (3.8)

Apparently if the ratio f/ρ is close to a constant, the variance is reduced.

It is tempting to take ρ = |f| and sample directly from the integrand. However in practice evaluations of the integrand are typically expensive. Therefore a better strategy is to build an approximate density in the product form, ρ(x, y, …, z) = ρ_{x}(x)ρ_{y}(y) … ρ_{z}(z), and then sample from this approximate density. A popular routine of this sort is called VEGAS.

3.4 Stratified sampling

Stratified sampling is a generalisation of the recursive adaptive integration algorithm to random quadratures in multi-dimensional spaces.

The ordinary “dividing by two” strategy does not work for multi-dimensional integrations as the number of sub-volumes grows way too fast to keep track

Table 3.2: Recursive stratified sampling algorithm

sample N random p o i n t s with p l a i n Monte Carlo ;
estimate the average and the e r r o r ;
IF the e r r o r i s a c c e p t a b l e :
RETURN the average and the e r r o r ;
ELSE :
FOR EACH dimension :
subdivide the volume in two along the dimension ;
estimate the sub−v a r i a n c e s in the two sub−volumes ;
pick the dimension with the l a r g e s t sub−v a r i a n c e ;
subdivide the volume in two along t h i s dimension ;
dispatch two r e c u r s i v e c a l l s to each o f the sub−volumes ;
estimate the grand average and grand e r r o r ;
RETURN the grand average and grand e r r o r ;

of. Instead one estimates along which dimension a subdivision should bring the most dividends and only subdivides along this dimension. Such strategy is called recursive stratified sampling. A simple variant of this algorithm is presented in Table 3.2.

Figure 3.1: Stratified sample of a discontinuous function

f(x, y) = {_{0 otherwise}^{1 if x^{2}+y^{2}<0.9^{2}}

built with the algorithm in Table 3.2.

In a stratified sample the points are concentrated in the regions where the variance of the function is largest, see an illustration in Figure 3.1.

Figure 3.2: Typical distributions of pseudo-random points (left), and quasirandom low-discrepancy points: lattice (center) and base-2/3 Halton (right) sequences. The first thousand points are plotted in each case.

3.5 Quasi-random (low-discrepancy) sampling

Pseudo-random sampling has high discrepancy^{1}: it typically creates regions with hight density of points and other regions with low density of points, see an illustration on Figure 3.2 (left). With pseudo-random sampling there is a finite probability that all the N points would fall into one half of the region and none into the other half.

Quasi-random sequences avoid this phenomenon by distributing points in a highly correlated manner with a specific requirement of low discrepancy, see Figure 3.2 for an example. Quasi-random sampling is like a computation on a grid where the grid constant must not be known in advance as the grid is ever gradually refined and the points are always distributed uniformly over the region. The computation can be stopped at any time.

By placing points more evenly than at random, the quasi-random sequences try to improve on the 1/\sqrt{N} convergence rate of pseudo-random sampling.

The central limit theorem does not apply in this case as the points are not statistically independent. Therefore the variance can not be used as an estimate of the error. The error estimation is actually not trivial. In practice one can employ two different sequences and use the difference in the resulting integrals as an error estimate.

3.5.1 Van der Corput and Halton sequences

A van der Corput sequence is a low-discrepancy sequence over the unit interval. It is constructed by reversing the base-b representation of the sequence of natural numbers (1, 2, 3, …). For example, the decimal van der Corput sequence begins as

0.1, 0.2, 0.3, …, 0.8, 0.9, 0.01, 0.11, 0.21, 0.31, …, 0.91, 0.02, 0.12, …. (3.9)

1discrepancy is a measure of how unevenly the points are distributed over the region.

In a base-b representation a natural number n with s digits {d_{i} | i = 1 … s, 0 ≤ d_{i} < b} is given as

n = ∑_{k=1}^{s} d_{k}b^{k−1} .  (3.10)

The corresponding base-b van der Corput number q_{b}(n) is then given as

q_{b}(n) = ∑_{k=1}^{s} d_{k}b^{−k} .  (3.11)

Here is a C implementation of this algorithm,

double corput ( int n , int base ){
double q=0, bk=(double )1/ base ;
while (n>0){ q += ( n % base )∗ bk ; n /= base ; bk /= base ; }
return q ; }

The van der Corput numbers of any base are uniformly distributed over the unit interval. They also form a dense set in the unit interval: there exists a subsequence of the van der Corput sequence which converges to any given real number in [0, 1].

The Halton sequence is a generalization of the van der Corput sequence to d-dimensional spaces. One chooses a set of coprime bases b_{1}, …, b_{d} and then for each dimension i generates a van der Corput sequence with its own base b_{i}. The n-th Halton d-dimentional point x in the unit volume is then given as

x_{b_{1},...,b_{d}}(n) = {q_{b_{1}}(n), …, q_{b_{d}}(n)} .  (3.12)

Here is a C implementation which calls the corput function listed above,

#include<a s s e r t . h>
void halton ( int n , int d , double ∗x ){
int base [ ] = { 2 , 3 , 5 , 7 , 1 1 , 1 3 , 1 7 , 1 9 , 2 3 , 2 9 , 3 1 , 3 7 , 4 1 , 4 3 , 4 7 , 5 3 , 5 9 , 6 1 , 6 7 };
int maxd=19; a s s e r t ( d <= maxd ) ;
for ( int i =0; i<d ; i++) x [ i ]= corput (n , base [ i ] ) ; }

3.5.2 Lattice rules

In the simplest incarnation a lattice rule can be defined as follows. Let the generating vector z = {α_{i} | i = 1, …, d}—where d is the dimension of the integration space—be a set of cleverly chosen irrational numbers. Then the n-th point (in the unit volume) of the sequence is given as

x(n) = frac(nz) ≡ {frac(nα_{1}), …, frac(nα_{d})} ,  (3.13)

where frac(x) is the fractional part of x.

An implementation of this algorithm in C is given in Table 3.3 and an illustration of such sequence is shown on Figure 3.2 (center).

3.6. IMPLEMENTATIONS 33

Table 3.3: Lattice low-disrepancy quasi-random sequence in C.

#include <math . h>
#include <a s s e r t . h>
#include <gc . h> // garbage c o l l e c t o r
#define f r a c l ( x ) ( ( x)− f l o o r l ( x ) )
#define r e a l long double
void l a t t i c e ( int d , double ∗x ){
static int dim=0, n=−1; static r e a l ∗ alpha ; int i ;
i f (d<0){ // d<0 i s the s i g n a l to ( re −) i n i t i a l i z e the l a t t i c e
dim=−d ; n=0; alpha=( r e a l ∗)GC MALLOC( dim∗ sizeof ( r e a l ) ) ;
for ( i =0; i<dim ; i++) alpha [ i ]= f r a c l ( s q r t l ( ( r e a l )13∗( i +1))); }
else {
a s s e r t ( d=dim && n>=0);
n++;for ( i =0; i<dim ; i++)x [ i ]= f r a c l ( n∗ alpha [ i ] ) ; }

}

3.6 Implementations

34 CHAPTER 3. MONTE CARLO INTEGRATION Table 3.4: Javascript implementation of the stratified sampling algorithm.

Array . prototype . i t e r a t o r =f u n c t i o n ( ) {
for ( var i =0; i<t h i s . length ; i++) y i e l d i ; }
f u n c t i o n s t r a t a ( f , a , b , acc , eps ,N, o l d s t a t ,V)

{

var randomx=f u n c t i o n ( a , b ) [ a [ i ]+Math . random ( ) ∗ ( b [ i ]−a [ i ] ) for ( i in a ) ]
var range =f u n c t i o n ( n ) { for ( var i =0; i<n ; i++) y i e l d i }
var s t a t s =f u n c t i o n ( xs ){
var n=xs . length ;
var average=xs . reduce ( f u n c t i o n ( a , b ) a+b , 0 ) / n ;
var v a r i a n c e=xs . reduce ( f u n c t i o n ( a , b ) a+b∗b , 0 ) / n−average ∗ average ;
return [ average , variance , n ] ;

}

i f ( typeof (N)==” undefined ” )N=42;
i f ( typeof ( o l d s t a t)==” undefined ” ){ // f i r s t c a l l : s e t t i n g up o l d s t a t s
var V=1; for ( l e t k in a ) V∗=(b [ k]−a [ k ] ) ;
var xs =[randomx ( a , b ) for ( i in range (N ) ) ] ;
var ys =[ f ( xs [ i ] ) for ( i in xs ) ] ;
var o l d s t a t=s t a t s ( ys ) ;

}

var xs =[randomx ( a , b ) for ( i in range (N) ) ] // new p o i n t s
var ys =[ f ( xs [ i ] ) for ( i in xs ) ]  // new f u n c t i o n v a l u e s
var [ average , variance ,]= s t a t s ( ys ) // average and variance
var [ oaverage , ovariance , oN]= o l d s t a t
var i n t e g=V∗( average ∗N+oaverage ∗oN) / (N+oN)
// i n t e g r a l and e r r o r
var e r r o r=V∗Math . s q r t ( ( v a r i a n c e ∗N+ovariance ∗oN) / (N+oN) / (N+oN) )
i f ( error <acc+eps ∗Math . abs ( i n t e g ) ) return [ integ , e r r o r ] ; // done
else { // not done : need to d i s p a t c h a r e c u r s i v e c a l l
var vmax=−1, kmax=0
for ( l e t k in a ){ // l o o k in a l l dimensions with i s b e s t to b i s e c t
var [ al , vl , nl ]=
s t a t s ( [ ys [ i ] for ( i in xs ) i f ( xs [ i ] [ k]< ( a [ k]+b [ k ] ) / 2 ) ] )
var [ ar , vr , nr ]=
s t a t s ( [ ys [ i ] for ( i in xs ) i f ( xs [ i ] [ k]>=(a [ k]+b [ k ] ) / 2 ) ] )
var v=Math . abs ( al−ar ) // t a k e the one with l a r g e s t v a r i a t i o n
i f ( v>vmax){ // remember the v a l u e s
vmax=v ; kmax=k ;
var o l d s t a t l =[ al , vl , nl ] , o l d s t a t r =[ar , vr , nr ]
}
} // now d i s p a t c h two r e c u r s i v e c a l l s
var a2=a . s l i c e ( ) ; a2 [ kmax]=( a [ kmax]+b [ kmax ] ) / 2
var b2=b . s l i c e ( ) ; b2 [ kmax]=( a [ kmax]+b [ kmax ] ) / 2
var [ i1 , e1 ]= s t r a t a ( f , a , b2 , acc / 1 . 4 1 4 , eps ,N, o l d s t a t l ,V/2)
var [ i2 , e2 ]= s t r a t a ( f , a2 , b , acc / 1 . 4 1 4 , eps ,N, o l d s t a t r ,V/2)
return [ i 1+i2 , Math . s q r t ( e1 ∗ e1+e2 ∗ e2 ) ] // return r e s u l t s

} }

Ordinary differential equations

4.1 Introduction

Ordinary differential equations (ODE) are generally defined as differential equations in one variable where the highest order derivative enters linearly. Such equations invariably arise in many different contexts throughout mathematics and science as soon as changes in the phenomena at hand are considered, usually with respect to variations of certain parameters.

Systems of ordinary differential equations can be generally reformulated as systems of first-order ordinary differential equations,

y′(x) = f(x, y) ,  (4.1)

where y′ ≐ dy/dx, and the variables y and the right-hand side function f(x, y) are understood as column-vectors. For example, a second order differential equation in the form

u′′ = g(x, u, u′)  (4.2)

can be rewritten as a system of two first-order equations,

\{ \begin{matrix}y′_{1} = y_{2} \\ y′_{2} = g(x, y_{1}, y_{2}) \end{matrix} ,  (4.3)

using the variable substitution y_{1} = u, y_{2} = u′.

In practice ODEs are usually supplemented with boundary conditions which pick out a certain class or a unique solution of the ODE. In the following we shall mostly consider initial value problems: ODE with the boundary condition in the form of an initial condition at a given point a,

y(a) = y_{0} .  (4.4)

The problem then is to find the value of the solution y at some other point b. Finding a solution to an ODE is often referred to as integrating the ODE.

An ODE integration algorighm typically advances the solution from the initial point a to the final point b in a number of discrete steps

{x_{0} ≐ a, x_{1}, …, x_{n−1}, x_{n} ≐ b}.  (4.5)

An efficient algorithm tries to integrate an ODE using as few steps as possible under the constraint of the given accuracy goal. For this purpose the algorthm should continuously adjust the step-size during the integration, using few larger steps in the regions where the solution is smooth and perhaps many smaller steps in more treacherous regions.

Typically, an adaptive step-size ODE integrator is implemented as two routines. One of them—called driver—monitors the local errors and tolerances and adjusts the step-sizes. To actually perform a step the driver calls a separate routine—the stepper—which advances the solution by one step, using one of the many available algorithms, and estimates the local error. The GNU Scientific Library, GSL, implements about a dozen of different steppers and a tunable adaptive driver.

In the following the chapter describes several of the popular driving algorithms and stepping methods for solving initial value ODE problems.

4.2 Error estimate

In an adaptive step-size algorithm the stepping routine must provide an estimate of the integration error, upon which the driver bases its strategy to determine the optimal step-size for a user-specified accuracy goal.

A stepping method is generally characterized by its order : a method has order p if it can integrate exactly an ODE where the solution is a polynomial of order p. In other words, for small h the error of the order-p method is O(h^{p+1}).

For sufficiently small steps the error δy of an integration step for a method of a given order p can be estimated by comparing the solution y_{full step} obtained with one full-step integration with a potentially more precise solution, y_{two half steps}, obtained with two consecutive half-step integrations,

δy = \frac{y_{full step} − y_{two half steps}}{2^{p+1} − 1} .  (4.6)

where p is the order of the algorithm used.

Indeed, if the step-size h is small, we can assume

δy_{full step} = Ch^{p+1} ,  (4.7)
δy_{two half steps} = C (\frac{h}{2})^{p+1} ,  (4.8)

where δy_{full step} and δy_{two half steps} are the errors of the full-step and two half-steps integrations, and C is an unknown constant. The two can be combined as

y_{full step} − y_{two half steps} = δy_{full step} − δy_{two half steps}
= C (\frac{h}{2})^{p+1} (2^{p+1} − 1) ,  (4.9)

from which it immedliately follows that

C (\frac{h}{2})^{p+1} = \frac{y_{full step} − y_{two half steps}}{2^{p+1} − 1} .  (4.10)

One has, of course, to take the potentially more precise y_{two half steps} as the solution y. Its error is then given as

δy ≐ δy_{two half steps} = C (\frac{h}{2})^{p+1} = \frac{y_{full step} − y_{two half steps}}{2^{p+1} − 1} , (4.11)

which had to be demonstrated. This prescription is often referred to as the Runge’s principle.

One drawback of the Runge’s principle is that the full-step and the two half-step calculations generally do not share evaluations of the right-hand side function f(x, y), and therefore many extra evaluations are needed to estimate the error.

An alternative prescription for error estimation is to make the same step-size integration using two methods of different orders, with the difference between the two solutions providing the estimate of the error. If the lower order method mostly uses the same evaluations of the right-hand side function—in which case it is called embedded in the higher order method— the error estimate does not need additional evaluations.

Predictor-corrector methods are naturally of embedded type: the correction — which generally increases the order of the method — itself can serve as the estimate of the error.

4.3 Runge-Kutta methods

Runge-Kutta methods are one-step methods which advance the solution over the current step using only the information gathered from withing the step itself. The solution y is advanced from the point x_{i} to x_{i+1} = x_{i} + h, where h is the step-size, using a one-step formula,

y_{i+1} = y_{i} + hk,  (4.12)

where y_{i+1} is the approximation to y(x_{i+1}), and the value k is chosen such that the method integrates exactly an ODE whose solution is a polynomial of the highest possible order.

The Runge-Kutta methods are distinguished by their order : a method has order p if it can integrate exactly an ODE where the solution is a polynomial of order p. In other words, for small h the error of the order-p method is O(h^{p+1}).

The first order Runge-Kutta method is the Euler’s method,

k = f(x_{0}, y_{0}) .  (4.13)

Second order Runge-Kutta methods advance the solution by an auxiliary evaluation of the derivative, e.g. the mid-point method,

k_{0} = f(x_{0}, y_{0}) ,
k_{1/2} = f(x_{0} + \frac{1}{2}h, y_{0} + \frac{1}{2}hk_{0}) ,
k = k_{1/2} ,  (4.14)

or the two-point method, also called the Heun’s method

k_{0} = f(x_{0}, y_{0}),
k_{1} = f(x_{0} + h, y_{0} + hk_{0}),
k = \frac{1}{2}(k_{0} + k_{1}) .  (4.15)

These two methods can be combined into a third order method,

k = \frac{1}{6}k_{0} + \frac{4}{6}k_{1/2} + \frac{1}{6}k_{1} .  (4.16)

The most commont is the fourth-order method, which is called RK4 or simply the Runge-Kutta method,

k_{0} = f(x_{0}, y_{0}) ,
k_{1} = f(x_{0} + \frac{1}{2}h, y_{0} + \frac{1}{2}hk_{0}) ,
k_{2} = f(x_{0} + \frac{1}{2}h, y_{0} + \frac{1}{2}hk_{1}) ,
k_{3} = f(x_{0} + h, y_{0} + hk_{2}) ,
k = \frac{1}{6}k_{0} + \frac{1}{3}k_{1} + \frac{1}{3}k_{2} + \frac{1}{6}k_{3} .  (4.17)

A general Runge-Kutta method can be written as

y_{n+1} = y_{n} + ∑_{i=1}^{s} b_{i}K_{i} ,  (4.18)

where

K_{1} = hf(x_{n}, y_{n}) ,
K_{2} = hf(x_{n} + c_{2}h, y_{n} + a_{21}K_{1}) ,
K_{3} = hf(x_{n} + c_{3}h, y_{n} + a_{31}K_{1} + a_{32}K_{2}) ,  (4.19)
⋮
K_{s} = hf(x_{n} + c_{s}h, y_{n} + a_{s1}K_{1} + a_{s2}K_{2} + ⋯ + a_{s,s−1}K_{s−1}) .

The upper case K_{i} are simply the lower case k_{i} multiplied by the step-size h.

To specify a particular Runge-Kutta method one needs to provide the coefficients {a_{ij}|1 ≤ j < i ≤ s}, {b_{i}|i = 1..s} and {c_{i}|i = 1..s}. The matrix [a_{ij}] is called the Runge-Kutta matrix, while the coefficients b_{i} and c_{i} are known as the weights and the nodes. These data are usually arranged in the so called Butcher’s tableau,

\begin{aligned}0 \\ c_{2} & a_{21} \\ c_{3} & a_{31} & a_{32} \\ ⋮ & ⋮ & & ... \\ c_{s} & a_{s1} & a_{s2} & ⋯ & a_{s,s−1} \\ & b_{1} & b_{2} & ⋯ & b_{s−1} & b_{s} \end{aligned} .  (4.20)

For example, the Butcher’s tableau for the RK4 method is

\begin{aligned}0 \\ 1/2 & 1/2 \\ 1/2 & 0 & 1/2 \\ 1 & 0 & 0 & 1 \\ & 1/6 & 1/3 & 1/3 & 1/6 \end{aligned} .  (4.21)

4.3.1 Embeded methods with error estimates

The embedded Runge-Kutta methods—in addition to advancing the solution by one step—produce an estimate of the local error of the step. This is done by having two methods in the tableau, one with a certain order p and another one with order p − 1. The difference bitween the two method gives the the estimate of the local error. The lower order method is embedded in the higher order method, that is, it uses the same K-values. This allows a very effective estimate of the error.

The embedded lower order method is written as

y_{n+1}^{∗} = y_{n} + ∑_{i=1}^{s} b_{i}^{∗}K_{i} ,  (4.22)

where K_{i} are the same as for the higher order method. The error estimate is then given as

e_{n} = y_{n+1} − y_{n+1}^{∗} = ∑_{i=1}^{s}(b_{i} − b_{i}^{∗})K_{i} .  (4.23)

The Butcher’s tableau for this kind of method is extended by one row to give the values of b_{i}^{∗}.

The simplest embedded methods are Heun-Euler method,

\begin{aligned}0 \\ 1 & 1 \\ & 1/2 & 1/2 \\ & 1 & 0 \end{aligned} ,  (4.24)

and midpoint-Euler method,

\begin{aligned}0 \\ 1/2 & 1/2 \\ & 0 & 1 \\ & 1 & 0 \end{aligned} ,  (4.25)

which both combine methods of orders 2 and 1. An implementation of the midpoint-Euler method in C is shown in Table 4.1.

Table 4.1: Embedded Runge-Kutta midpoint-Euler stepper with error estimate.

void rkstep12 ( void f ( int n , double x , double∗ y , double∗ dydx ) ,
int n , double x , double∗ y , double h , double∗ yout , double∗ e r r )
{ double k1 [ n ] , k2 [ n ] , ytmp [ n ] ; //VLA: −s t d=c99
f (n , x , y , k1 ) ;
for ( int i = 0 ; i<n ; i++) ytmp [ i ]=y [ i ]+1./2∗ h∗k1 [ i ] ;
f (n , x+1./2∗h , ytmp , k2 ) ;
for ( int i = 0 ; i<n ; i ++){ yout [ i ]=y [ i ]+h∗k2 [ i ] ;
e r r [ i ]=h∗( k2 [ i ]−k1 [ i ] ) / 2 ; } // o p t i m i s t i c

}

The Bogacki-Shampine method [2] combines methods of orders 3 and 2,

\begin{aligned}0 \\ 1/2 & 1/2 \\ 3/4 & 0 & 3/4 \\ 1 & 2/9 & 1/3 & 4/9 \\ & 2/9 & 1/3 & 4/9 & 0 \\ & 7/24 & 1/4 & 1/3 & 1/8 \end{aligned} .  (4.26)

Bogacki and Shampine argue that their method has better stability properties and actually outperforms higher order methods at lower accuracy goal calculations. This method has the FSAL—first same as last—property: the value k_{4} at one step equals k_{1} at the next step; thus only three function evaluations are needed per step. Table 4.2 shows a simple implementation which does not utilise this property for the sake of presentational clarity.

Table 4.2: Embedded Runge-Kutta Bogacki-Shampine stepper

void rkstep23 ( void f ( int n , double x , double∗ y , double∗ dydx ) ,
int n , double x , double∗ y , double h , double∗ yout , double∗ e r r )
{ double k1 [ n ] , k2 [ n ] , k3 [ n ] , k4 [ n ] , ytmp [ n ] ; //VLA: −s t d=c99
f (n , x , y , k1 ) ;
for ( int i =0; i<n ; i++) ytmp [ i ]=y [ i ]+1./2∗ k1 [ i ] ∗ h ;
f (n , x+1./2∗h , ytmp , k2 ) ;
for ( int i =0; i<n ; i++) ytmp [ i ]=y [ i ]+3./4∗ k2 [ i ] ∗ h ;
f (n , x+3./4∗h , ytmp , k3 ) ;
for ( int i =0; i<n ; i++)
ytmp [ i ]=y [ i ]+(2./9∗ k1 [ i ]+1./3∗ k2 [ i ]+4./9∗ k3 [ i ] ) ∗ h ;
f (n , x+h , ytmp , k4 ) ;
for ( int i = 0 ; i<n ; i++)
{ yout [ i ]=y [ i ]+(2./9∗ k1 [ i ]+1./3∗ k2 [ i ]+4./9∗ k3 [ i ] ) ∗ h ;
ytmp [ i ]=y [ i ]+(7./24∗ k1 [ i ]+1./4∗ k2 [ i ]+1./3∗ k3 [ i ]+1./8∗ k4 [ i ] ) ∗ h ;
e r r [ i ]= yout [ i ]−ytmp [ i ] ;

} }

The Runge-Kutta-Fehlberg method [5] — called RKF45 — implemented in the renowned rkf45 Fortran routine, has two methods of orders 5 and 4,

0

1/4  1/4
3/8 3/32 9/32
12/13 1932/2197 −7200/2197 7296/2197
1  439/216  −8  3680/513 −845/4104
1/2  −8/27  2  −3544/2565 1859/4104 −11/40
16/135  0  6656/12825 28561/56430 −9/50 2/55
25/216  0  1408/2565 2197/4104 −1/5  0

4.4 Multistep methods

Multistep methods try to use the information about the function gathered at the previous steps. They are generally not self-starting as there are no previous points at the start of the integration. The first step must be done with a one-step method like Runge-Kutta.

A number of multistep have been devised (and named after different mathematicians); we shall only consider a few simple ones here to get the idea of how it works.

4.4.1 Two-step method

Given the previous point, (x_{i−1}, y_{i−1}), in addition to the current point (x_{i}, y_{i}), the sought function y can be approximated in the vicinity of the point x_{i} as

ȳ(x) = y_{i} + y′_{i} · (x − x_{i}) + c̄ · (x − x_{i})^{2},  (4.27)

where y′_{i} = f(x_{i}, y_{i}) and the coefficient c̄ can be found from the condition ȳ(x_{i−1}) = y_{i−1}, which gives

c̄ = \frac{y_{i−1} − y_{i} + y′_{i} · (x_{i} − x_{i−1})}{(x_{i} − x_{i−1})^{2}} .  (4.28)

The value of the function at the next point, x_{i+1} ≐ x_{i} + h, can now be estimated as ȳ(x_{i+1}) from (4.27).

The error of this second-order two-step stepper can be estimated by a comparison with the first-order Euler’s step, which is given by the linear part of (4.27). The correction term c̄h^{2} can serve as the error estimate,

δy = c̄h^{2} .  (4.29)

4.4.2 Two-step method with extra evaluation

One can further increase the order of the approximation (4.27) by adding a third order term,

¯ȳ(x) = ȳ(x) + ¯ d̄ · (x − x_{i})^{2}(x − x_{i−1}) .  (4.30)

The coefficient ¯ d̄ can be found from the matching condition at a certain point inside the invervarl, for example at half-step,

¯ȳ′(x_{1/2}) = f(x_{1/2}, ȳ(x_{1/2})) ≐ f̄_{1/2} ,  (4.31)

where x_{1/2} ≐ x_{i} + h/2. This gives

¯d̄ = \frac{f̄_{1/2} − y′_{i} − 2c̄ · (x_{1/2} − x_{i})}{2(x_{1/2} − x_{i})(x_{1/2} − x_{i−1}) + (x_{1/2} − x_{i})^{2}} .  (4.32)

The error estimate at the point x_{i+1} ≐ x_{0} +h is again given as the difference between the higher and the lower order methods,

δy = ¯ ȳ(x_{i+1}) − ȳ(x_{i+1}) .  (4.33)

4.5 Predictor-corrector methods

A predictor-corrector method uses extra iterations to improve the solution. It is an algorithm that proceeds in two steps. First, the predictor step calculates a rough approximation of y(x + h). Second, the corrector step refines the initial approximation. Aditionally the corrector step can be repeated in the hope that this achieves an even better approximation to the true solution.

For example, the two-point Runge-Kutta method (4.15) is as actually a predictor-corrector method, as it first calculates the prediction ỹ_{i+1} for y(x_{i+1}),

ỹ_{i+1} = y_{i} + hf(x_{i}, y_{i}) ,  (4.34)

and then uses this prediction in a correction step,

ˇỹ_{i+1} = y_{i} + h\frac{1}{2} (f(x_{i}, y_{i}) + f(x_{i+1}, ỹ_{i+1})) .  (4.35)

4.5.1 Two-step method with correction

Similarly, one can use the two-step approximation (4.27) as a predictor, and then improve it by one order with a correction step, namely

ˇȳ(x) = ȳ(x) + ď · (x − x_{i})^{2}(x − x_{i−1}).  (4.36)

The coefficient ď can be found from the condition ˇ ȳ′(x_{i+1}) = f̄_{i+1}, where f̄_{i+1} ≐ f(x_{i+1}, ȳ(x_{i+1})),

ď = \frac{f̄_{i+1} − y′_{i} − 2c̄ · (x_{i+1} − x_{i})}{2(x_{i+1} − x_{i})(x_{i+1} − x_{i−1}) + (x_{i+1} − x_{i})^{2}} .  (4.37)

Equation (4.36) gives a better estimate, y_{i+1} = ˇ ȳ(x_{i+1}), of the sought function at the point x_{i+1}. In this context the formula (4.27) serves as predictor, and (4.36) as corrector. The difference between the two gives an estimate of the error.

This method is equivalent to the two-step method with an extra evaluation where the extra evaluation is done at the full step.

4.6 Adaptive step-size control

Let tolerance τ be the maximal accepted error consistent with the required accuracy to be achieved in the integration of an ODE. Suppose the inegration is done in n steps of size h_{i} such that ∑_{i=1}^{n} h_{i} = b−a. Under assumption that the errors at the integration steps are random and statistically uncorrelated, the local tolerance τ_{i} for the step i has to scale as the square root of the step-size,

τ_{i} = τ\sqrt{\frac{h_{i}}{b − a}}.  (4.38)

Indeed, if the local error e_{i} on the step i is less than the local tolerance, e_{i} ≤ τ_{i}, the total error E will be consistent with the total tolerance τ,

E ≈ \sqrt{∑_{i=1}^{n} e_{i}^{2}} ≤ \sqrt{∑_{i=1}^{n} τ_{i}^{2}} = τ\sqrt{∑_{i=1}^{n} \frac{h_{i}}{b − a}} = τ .  (4.39)

The current step h_{i} is accepted if the local error e_{i} is smaller than the local tolerance τ_{i}, after which the next step is attempted with the step-size adjusted according to the following empirical prescription [4],

h_{i+1} = h_{i} × (\frac{τ_{i}}{e_{i}} )^{Power} × Safety,  (4.40)

where Power ≈ 0.25 and Safety ≈ 0.95.

If the local error is larger than the local tolerance the step is rejected and a new step is attempted with the step-size adjusted according to the same prescription (4.40).

One simple prescription for the local tolerance τ_{i} and the local error e_{i} to be used in (4.40) is

τ_{i} = (ϵ∥y_{i}∥ + δ)\sqrt{\frac{h_{i}}{b − a}} , e_{i} = ∥δy_{i}∥ ,  (4.41)

where δ and ϵ are the required absolute and relative precision and δy_{i} is the estimate of the integration error at the step i.

A more elaborate prescription considers components of the solution separately,

(τ_{i})_{k} = (ϵ|(y_{i})_{k}| + δ)\sqrt{\frac{h_{i}}{b − a}} , (e_{i})_{k} = |(δy_{i})_{k}| ,  (4.42)

where the index k runs over the components of the solution. In this case the step acceptence criterion also becomes component-wise: the step is accepted, if

∀k : (e_{i})_{k} < (τ_{i})_{k} .  (4.43)

The factor τ_{i}/e_{i} in the step adjustment formula (4.40) is then replaced by

\frac{τ_{i}}{e_{i}} → min_{k} \frac{(τ_{i})_{k}}{(e_{i})_{k}} .  (4.44)

Yet another refinement is to include the derivatives y′ of the solution into the local tolerance estimate, either overally,

τ_{i} = (ϵα∥y_{i}∥ + ϵβ∥y′_{i}∥ + δ)\sqrt{\frac{h_{i}}{b − a}} ,  (4.45)

or commponent-wise,

(τ_{i})_{k} = (ϵα|(y_{i})_{k}| + ϵβ|(y′_{i})_{k}| + δ)\sqrt{\frac{h_{i}}{b − a}} .  (4.46)

The weights α and β are chosen by the user.

Table (4.3) shows a simple implementation of the described algorithm.

4.6. ADAPTIVE STEP-SIZE CONTROL 45

Table 4.3: An ODE driver with adaptive step-size control in C.

#include <math . h>
#include <s t d l i b . h>
#include <s t d i o . h>
int d r i v e (
void f ( int n , double x , double∗ y , double∗ dydx ) , int n , double∗ x l i s t ,
double∗∗ y l i s t , double b , double h , double acc , double eps , int max ,
void step ( void f ( int n , double x , double∗ y , double∗ dydx ) ,
int n , double x , double∗ y , double h , double∗ yout , double∗ dy )

) {

int i =0; double a=x l i s t [ 0 ] ; // s t a r t i n g p o i n t
double dy [ n ] , y1 [ n ] ; // w i l l s t o r e d e l t a y and y { i +1}
while ( x l i s t [ i ]<b ) // do u n t i l b i s reached

{

double x=x l i s t [ i ] , ∗y=y l i s t [ i ] ; // current x and y
i f ( x+h>b ) h=b−x ; // make sure we land on b
step ( f , n , x , y , h , y1 , dy ) ; // the s t e p
double sum=0; for ( int k=0;k<n ; k++) sum+=dy [ k ] ∗ dy [ k ] ;
double e r r=s q r t (sum ) ; // l o c a l e r r o r
sum=0; for ( int k=0;k<n ; k++) sum+=y1 [ k ] ∗ y1 [ k ] ;
double normy=s q r t (sum ) ;
double t o l =(normy∗ eps+acc )∗ s q r t ( h/(b−a ) ) ; // l o c a l t o l e r a n c e :
i f ( t o l >e r r ) // then accept s t e p and s t o r e x and y

{

i ++; i f ( i=max) return −i ; // s t o r a g e f i l l e d , a b o r t : (
x l i s t [ i ]=x+h ; for ( int k=0;k<n ; k++) y l i s t [ i ] [ k]=y1 [ k ] ;
} // a d j u s t the s t e p :
i f ( err >0) h=fmin (2∗h , h∗pow( t o l / err , 0 . 2 5 ) ∗ 0 . 9 5 ) ; else h = 2∗h ;

}

return i +1;

}

46 CHAPTER 4. ORDINARY DIFFERENTIAL EQUATIONS

Systems of linear equations

5.1 Introduction

A system of linear equations (or linear system) is a collection of linear equations involving the same set of unknown variables. A general system of n linear equations with m unknowns can be written as

\{ \begin{matrix}A_{11}x_{1} & + & A_{12}x_{2} & + & ⋯ & + & A_{1m}x_{m} & = & b_{1} \\ A_{21}x_{1} & + & A_{22}x_{2} & + & ⋯ & + & A_{2m}x_{m} & = & b_{2} \\ ⋮ & & ⋮ & & & & ⋮ & & ⋮ \\ A_{n1}x_{1} & + & A_{n2}x_{2} & + & ⋯ & + & A_{nm}x_{m} & = & b_{n} \end{matrix} ,  (5.1)

where x_{1}, x_{2}, …, x_{m} are the unknown variables, A_{11}, A_{12}, …, A_{nm} are the (constant) coefficients of the system, and b_{1}, b_{2}, …, b_{n} are the (constant) right-hand side terms.

The system can be equivalently written in the matrix form,

Ax = b ,  (5.2)

where A ≐ {A_{ij}} is the n × m matrix of the coefficients, x ≐ {x_{j}} is the size-m column-vector of the unknown variables, and b ≐ {b_{i}} is the size-n column-vector of right-hand side terms.

A solution to a linear system is a set of values for the variables x which satisfies all equations.

Systems of linear equations occur quite regularly in applied mathematics. Therefore computational algorithms for finding solutions of linear systems are an important part of numerical methods.

A system of non-linear equations can often be approximated by a linear system – a helpful technique (called linearization) in creating a mathematical model of an otherwise a more complex system.

If m = n, the matrix A is called square. A square system has a unique solution if A is invertible.

5.2 Triangular systems

An efficient algorithm to solve a square system of linear equations numerically is to transform the original system into an equivalent triangular system,

Ty = c ,  (5.3)

where T is a triangular matrix – a special kind of square matrix where the matrix elements either below (lower triangular) or above (upper triangular) the main diagonal are zero.

Indeed, an upper triangular system Uy = c can be easily solved by back-substitution,

y_{i} = \frac{1}{U_{ii}} (c_{i} − ∑_{k=i+1}^{n} U_{ik}y_{k}) , i = n, n − 1, …, 1 ,  (5.4)

where one first computes y_{n} = b_{n}/U_{nn}, then substitutes back into the previous equation to solve for y_{n−1}, and repeats through y_{1}.

Here is a C-function implementing in-place^{1} back-substitution^{2}:

#include<matrix . h>
#include<a s s e r t . h>
void backsub ( matrix ∗U, v e c t o r ∗b){
a s s e r t (U−>s i z e 1=U−>s i z e 2 && U−>s i z e 2=b−>s i z e ) ;
for ( int i=b−>s i z e −1; i >=0;i −−){
double s=v e c t o r g e t (b , i ) ;
for ( int k=i +1;k<b−>s i z e ; k++) s−=matrix get (U, i , k )∗ v e c t o r g e t (b , k ) ;
v e c t o r s e t (b , i , s / matrix get (U, i , i ) ) ;

} }

For a lower triangular system Ly = c the equivalent procedure is called forward-substitution,

y_{i} = \frac{1}{L_{ii}} (c_{i} − ∑_{k=1}^{i−1} U_{ik}y_{k}) , i = 1, 2, …, n .  (5.5)

5.3 Reduction to triangular form

Popular algorithms for reducing a square system of linear equations to a triangular form are LU-decomposition and QR-decomposition.

1here in-place means the right-hand side c is replaced by the solution y.

2the functions vector get, vector set, and matrix get are assumed to implement fetching and setting the vector and matrix elements.

5.3.1 QR-decomposition

QR-decomposition is a factorization of a matrix into a product of an orthogonal matrix Q, such that Q^{⊺}Q = 1, where ^{⊺} denotes transposition, and a right triangular matrix R,

A = QR .  (5.6)

QR-decomposition can be used to convert a linear system Ax = b into the triangular form (by multiplying with Q^{⊺} from the left),

Rx = Q^{⊺}b ,  (5.7)

which can be solved directly by back-substitution.

QR-decomposition can also be performed on non-square matrices with few long columns. Generally speaking a rectangular n × m matrix A can be represented as a product, A = QR, of an orthogonal n × m matrix Q, Q^{⊺}Q = 1, and a right-triangular m × m matrix R.

QR-decomposition of a matrix can be computed using several methods, such as Gram-Schmidt orthogonalization, Householder transformation [10], or Givens rotation [6].

Gram-Schmidt orthogonalization

Gram-Schmidt orthogonalization is an algorithm for orthogonalization of a set of vectors in a given inner product space. It takes a linearly independent set of vectors A = {a_{1}, …, a_{m}} and generates an orthogonal set Q = {q_{1}, …, q_{m}} which spans the same subspace as A. The algorithm is given as

for i = 1 to m :

q_{i} ← a_{i}/∥a_{i}∥

for j = i + 1 to m : a_{j} ← a_{j} − ⟨q_{i}|a_{j}⟩ q_{i}

where ⟨a|b⟩ is the inner product of two vectors, and ∥a∥ ≐ \sqrt{⟨a|a⟩} is the vector’s norm. This variant of the algorithm, where all remaining vectors a_{j} are made orthogonal to q_{i} as soon as the latter is calculated, is considered to be numerically stable and is referred to as stabilized or modified.

Stabilized Gram-Schmidt orthogonalization can be used to compute QRdecomposition of a matrix A by orthogonalization of its column-vectors a_{i} with the inner product

⟨a|b⟩ = a^{⊺}b ≡ ∑_{k=1}^{n}(a)_{k}(b)_{k} ,  (5.8)

where n is the length of column-vectors a and b, and (a)_{k} is the kth element of the column-vector, Table 5.1: QR-decomposition in C using stabilized Gram-Schmidt orthogonalization.

#include <matrix . h>
#include <math . h>
#include <s t d i o . h>
void qrdec ( matrix ∗ A, matrix ∗ R){ // A i s r e p l a c e d with Q, R i s f i l l e d
for ( int i =0; i<A−>s i z e 2 ; i ++){
v e c t o r ∗ a i = matrix get column (A, i ) ;
double r = s q r t ( v e c t o r d o t ( ai , a i ) ) ;
m a t r i x s e t (R, i , i , r ) ;
v e c t o r s c a l e ( ai , 1 / r ) ;
for ( int j=i +1; j<A−>s i z e 2 ; j ++){
v e c t o r ∗ a j=matrix get column (A, j ) ;
double s=v e c t o r d o t ( ai , a j ) ;
vector add ( aj ,−s , a i ) ;
m a t r i x s e t (R, i , j , s ) ; } } }

for i = 1 to m :

R_{ii} = \sqrt{a_{i}^{⊺}a_{i}} ; q_{i} = a_{i}/R_{ii}

  for j = i + 1 to m :

R_{ij} = q_{i}^{⊺}a_{j} ; a_{j} = a_{j} − q_{i}R_{ij} .

After orthogonalization the matrices Q = {q_{1} … q_{m}} and R are the sought orthogonal and right-triangular factors of matrix A.

The factorization is unique under requirement that the diagonal elements of R are positive. For a n × m matrix the complexity of the algorithm is O(m^{2}n).

Table (5.1) shows a C implementation of the stabilized Gram-Schmidt orthogonalization. The corresponding function to solve the equation

QRx = b  (5.9)

are shown in Table (5.2).

Householder transformation

A square matrix H of the form

H = 1 − \frac{2}{u^{⊺}u}uu^{⊺}  (5.10)

is called Householder matrix, where the vector u is called a Householder vector. Householder matrices are symmetric and orthogonal,

H^{⊺} = H , H^{⊺}H = 1 .  (5.11)

The transformation induced by the Householder matrix on a given vector a,

a → Ha ,  (5.12)

Table 5.2: C-functions to perform QR-backsubstitution of equation (5.9)

#include<matrix . h> // garbage c o l l e c t e d m a t r i x l i b r a r y
#include<a s s e r t . h>
void backsub ( matrix ∗, v e c t o r ∗) ;
v e c t o r ∗ qrback ( matrix ∗Q, matrix ∗R, v e c t o r ∗b){
v e c t o r ∗ x = matrixT times vector (Q, b ) ;
backsub (R, x ) ;
return x ;

}

void q r b a c k i n p l a c e ( matrix ∗Q, matrix ∗R, v e c t o r ∗b){
a s s e r t (Q−>s i z e 1=Q−>s i z e 2 ) ;
v e c t o r ∗ x = qrback (Q,R, b ) ;
for ( int i =0; i<b−>s i z e ; i++)v e c t o r s e t (b , i , v e c t o r g e t ( x , i ) ) ;

}

is called a Householder transformation or Householder reflection. The transformation changes the sign of the affected vector’s component in the u direction, or, in other words, makes a reflection of the vector about the hyperplane perpendicular to u, hence the name.

Householder transformation can be used to zero selected components of a given vector a. For example, one can zero all components but the first one, such that

Ha = γe_{1} ,  (5.13)

where γ is a number and e_{1} is the unit vector in the first direction. The factor γ can be easily calculated,

∥a∥^{2} ≐ a^{⊺}a = a^{⊺}H^{⊺}Ha = (γe_{1})^{⊺}(γe_{1}) = γ^{2} ,  (5.14)
⇒ γ = ±∥a∥ .  (5.15)

To find the Householder vector, we notice that

a = H^{⊺}Ha = H^{⊺}γe_{1} = γe_{1} − \frac{2(u)_{1}}{u^{⊺}u} u ,  (5.16)

⇒ \frac{2(u)_{1}}{u^{⊺}u} u = γe^{1} − a ,  (5.17)

where (u)_{1} is the first component of the vector u. One usually chooses (u)_{1} = 1 (for the sake of the possibility to store the other components of the Householder vector in the zeroed elements of the vector a) and stores the factor

\frac{2}{u^{⊺}u} ≡ τ  (5.18)

separately. With this convention one readily finds τ from the first component of equation (5.17),

τ = γ − (a)_{1} .  (5.19)

where (a)_{1} is the first element of the vector a. For the sake of numerical stability the sign of γ has to be chosen opposite to the sign of (a)_{1},

γ = −sign ((a)_{1}) ∥a∥ .  (5.20)

Finally, the Householder reflection, which zeroes all component of a vector a but the first, is given as

H = 1 − τuu^{⊺} , τ = −sign((a)_{1})∥a∥ − (a)_{1} , (u)_{1} = 1 , (u)_{i>1} = −\frac{1}{τ} (a)_{i} .

(5.21)

Now, a QR-decomposition of an n × n matrix A by Householder transformations can be performed in the following way:

1. Build the size-n Householder vector u_{1} which zeroes the sub-diagonal elements of the first column of matrix A, such that

H_{1}A = [ \begin{matrix}⋆ & ⋆ & … & ⋆ \\ 0 \\ ⋮ & & A_{1} \\ 0 \end{matrix} ] ,  (5.22)

where H_{1} = 1 − τ_{1}u_{1}u_{1}^{⊺} and where ⋆ denotes (generally) non-zero matrix elements. In practice one does not build the matrix H_{1} explicitly, but rather calculates the matrix H_{1}A in-place, consecutively applying the Householder reflection to columns the matrix A, thus avoiding computationally expensive matrix-matrix operations. The zeroed subdiagonal elements of the first column of the matrix A can be used to store the elements of the Householder vector u_{1} while the factor τ_{1} has to be stored separately in a special array. This is the storage scheme used by LAPACK and GSL.

2. Similarly, build the size-(n−1) Householder vector u_{2} which zeroes the sub-diagonal elements of the first column of matrix A_{1} from eq. (5.22). With the transformation matrix H_{2} defined as

H_{2} = [ \begin{matrix}1 & 0 & ⋯ & 0 \\ 0 \\ ⋮ & & 1 − τ_{2}u_{2}u_{2}^{⊺} \\ 0 \end{matrix} ] .  (5.23)

the two transformations together zero the sub-diagonal elements of the two first columns of matrix A,

H_{2}H_{1}A = [ \begin{matrix}⋆ & ⋆ & ⋆ & ⋯ & ⋆ \\ 0 & ⋆ & ⋆ & ⋯ & ⋆ \\ 0 & 0 \\ ⋮ & ⋮ & & A_{3} \\ 0 & 0 \end{matrix} ] ,  (5.24)

3. Repeating the process zero the sub-diagonal elements of the remaining columns. For column k the corresponding Householder matrix is

H_{k} = [ \begin{matrix}I_{k−1} & 0 \\ 0 & 1 − τ_{k}u_{k}u_{k}^{⊺} \end{matrix} ] ,  (5.25)

where I_{k−1} is a unity matrix of size k − 1, u_{k} is the size-(n-k+1) Householder vector that zeroes the sub-diagonal elements of matrix A_{k−1} from the previous step. The corresponding transformation step is

H_{k} … H_{2}H_{1}A = [ \begin{matrix}R_{k} & ⋆ \\ 0 & A_{k} \end{matrix} ] ,  (5.26)

where R_{k} is a size-k right-triangular matrix.

After n − 1 steps the matrix A will be transformed into a right triangular matrix,

H_{n−1} ⋯ H_{2}H_{1}A = R .  (5.27)

4. Finally, introducing an orthogonal matrix Q = H_{1}^{⊺}H_{2}^{⊺} … H_{n−1}^{⊺} and multiplying eq. (5.27) by Q from the left, we get the sought QRdecomposition,

A = QR .  (5.28)

In practice one does not explicitly builds the Q matrix but rather applies the successive Householder reflections stored during the decomposition.

Givens rotations

A Givens rotation is a transformation in the form

A → G(p, q, θ)A ,  (5.29)

where A is the transformed object—matrix of vector—and G(p, q, θ) is the Givens rotation matrix (also known as Jacobi rotation matrix): an orthogonal matrix in the form

G(p, q, θ) = [\begin{matrix}1 & ⋯ & 0 & ⋯ & 0 & ⋯ & 0 \\ ⋮ & ... & ⋮ & & ⋮ & & ⋮ \\ 0 & ⋯ & cos θ & ⋯ & sin θ & ⋯ & 0 \\ ⋮ & & ⋮ & ... & ⋮ & & ⋮ \\ 0 & ⋯ & − sin θ & ⋯ & cos θ & ⋯ & 0 \\ ⋮ & & ⋮ & & ⋮ & ... & ⋮ \\ 0 & ⋯ & 0 & ⋯ & 0 & ⋯ & 1 \end{matrix}] \begin{matrix}← row p \\ ← row q \end{matrix} . (5.30)

When a Givens rotation matrix G(p, q, θ) multiplies a vector, only elements with indices p and q are affected. Restricting our attention to these two elements, say a and b, the Givens rotation is given explicitely as

[\begin{matrix}cos θ & sin θ \\ − sin θ & cos θ \end{matrix}] [\begin{matrix}a \\ b \end{matrix}] = [\begin{matrix}a cos θ + b sin θ \\ −a sin θ + b cos θ \end{matrix}] .  (5.31)

Apparently, the rotation can zero the lowest element, if the angle θ is chosen as

tan θ = \frac{b}{a} ⇒ θ = atan2(b, a) .  (5.32)

A sequence of Givens rotations can zero all elements of a matrix below the main diagonal. Which obviously amounts to QR-decomposition of the matrix with Q being the product of all the applied Givens matrices.

Since each Givens rotation only affects two rows of the matrix it is possible to apply a set of rotations in parallel. Givens rotations are also more efficient on sparse matrices.

5.3.2 LU-decomposition

LU-decomposition is a factorization of a square matrix A into a product of a lower triangular matrix L and an upper triangular matrix U,

A = LU .  (5.33)

The linear system Ax = b after LU-decomposition of the matrix A becomes LUx = b and can be solved by first solving Ly = b for y and then Ux = y for x with two runs of forward and backward substitutions.

If A is an n × n matrix, the condition (5.33) is a set of n^{2} equations,

∑_{k=1}^{n} L_{ik}U_{kj} = A_{ij} |_{i,j=1...n} ,  (5.34)

for n^{2} + n unknown elements of the triangular matrices L and U. The decomposition is thus not unique.

Usually the decomposition is made unique by providing extra n conditions e.g. by the requirement that the elements of the main diagonal of the matrix L are equal one,

L_{ii} = 1 , i = 1 … n .  (5.35)

The system (5.34) with the extra conditions (5.35) can then be easily solved row after row using the Doolittle’s algorithm,

for i = 1 … n :

L_{ii} = 1

for j = i … n : U_{ij} = A_{ij} − ∑_{k<i} L_{ik}U_{kj}

for j = i + 1 … n : L_{ji} = \frac{1}{U_{ii}} (A_{ji} − ∑_{k<j} L_{jk}U_{ki})

In a slightly different Crout’s algorithm it is the matrix U that has unit diagonal elements,

for i = 1 … n :

U_{ii} = 1

for j = i … n : L_{ji} = A_{ji} − ∑_{k<i} L_{jk}U_{ki}

for j = i + 1 … n : U_{ij} = \frac{1}{L_{ii}} (A_{ji} − ∑_{k<j} L_{jk}U_{ki})

Without a proper ordering (permutations) in the matrix, the factorization may fail. For example, it is easy to verify that A_{11} = L_{11}U_{11}. If A_{11} = 0, then at least one of L_{11} and U_{11} has to be zero, which implies either L or U is singular, which is impossible if A is non-singular. This is however only a procedural problem. It can be removed by simply reordering the rows of A so that the first element of the permuted matrix is nonzero (or, even better, the largest in absolute value among all elements of the column below the diagonal). The same problem in subsequent factorization steps can be removed in a similar way. Such algorithm is refered to as partial pivoting. It requires an extra integer array to keep track of row permutations.

5.4 Determinant of a matrix

LU- and QR-decompositions allow O(n^{3}) calculation of the determinant of a square matrix. Indeed, for the LU-decomposition,

det A = det LU = det L det U = det U = ∏_{i=1}^{n} U_{ii} .  (5.36)

For the QR-decomposition

det A = det QR = det Q det R .  (5.37)

Since Q is an orthogonal matrix (det Q)^{2} = 1 and therefore

| det A| = | det R| = |∏_{i=1}^{n} R_{ii}| .  (5.38)

5.5 Matrix inverse

The inverse A^{−1} of a square n × n matrix A can be calculated by solving n linear equations

Ax_{i} = e_{i} |_{i=1...n} ,  (5.39)

where e_{i} is a column where all elements are equal zero except for the element number i, which is equal one; that is, columns e_{i} form a unity matrix. The matrix made of columns x_{i} is apparently the inverse of A.

A C implementation of this algorithm using QR-decomposition is shown in Table (5.3).

56 CHAPTER 5. SYSTEMS OF LINEAR EQUATIONS

Table 5.3: Matrix inverse with QR-decomposition in C

#include <matrix . h>
void q r b a c k i n p l a c e ( matrix ∗, matrix ∗, v e c t o r ∗) ;
void q r i n v e r s e ( matrix ∗Q, matrix ∗R, matrix ∗Ainv ){
m a t r i x s e t i d e n t i t y ( Ainv ) ;
for ( int i =0; i<Q−>s i z e 1 ; i++) q r b a c k i n p l a c e (Q,R, matrix get column ( Ainv , i ) ) ;

}

Ordinary least squares

6.1 Introduction

A system of linear equations is considered overdetermined if there are more equations than unknown variables. If all equations of an overdetermined system are linearly independent, the system has no exact solution.

A linear least-squares problem is the problem of finding an approximate solution to an overdetermined linear system. It often arises in applications where a theoretical model is fitted to experimental data.

6.2 Linear least-squares problem

Consider a linear system

Ac = b ,  (6.1)

where A is an n × m matrix, c is an m-component vector of unknown variables and b is an n-component vector of the right-hand side terms. If the number of equations n is larger than the number of unknowns m, the system is overdetermined and generally has no solution.

However, it is still possible to find an approximate solution—the one where Ac is only approximately equal b—in the sence that the Euclidean norm of the difference between Ac and b is minimized,

c : min_{c} ∥Ac − b∥^{2} .  (6.2)

The problem (6.2) is called the ordinary least-squares problem and the vector c that minimizes ∥Ac − b∥^{2} is called the least-squares solution.

6.3 Solution via QR-decomposition

The linear least-squares problem can be solved by QR-decomposition. The matrix A is factorized as A = QR, where Q is n×m matrix with orthogonal columns, Q^{T} Q = 1, and R is an m × m upper triangular matrix. The Euclidean norm ∥Ac − b∥^{2} can then be rewritten as

∥Ac − b∥^{2} = ∥QRc − b∥^{2}
= ∥Rc − Q^{T} b∥^{2} + ∥(1 − QQ^{T})b∥^{2} ≥ ∥(1 − QQ^{T})b∥^{2} .  (6.3)

The term ∥(1 − QQ^{T})b∥^{2} is independent of the variables c and can not be reduced by their variations. However, the term ∥Rc−Q^{T} b∥^{2} can be reduced down to zero by solving the m × m system of linear equations

Rc = Q^{T} b .  (6.4)

The system is right-triangular and can be readily solved by back-substitution. Thus the solution to the ordinary least-squares problem (6.2) is given by the solution of the triangular system (6.4).

6.4 Ordinary least-squares curve fitting

Ordinary least-squares curve fitting is a problem of fitting n (experimental) data points {x_{i}, y_{i} ± Δy_{i}}, where Δy_{i} are experimental errors, by a linear combination, F, of m functions {f_{k}(x) | k = 1, …, m},

F_{c}(x) = ∑_{k=1}^{m} c_{k}f_{k}(x) .  (6.5)

The objective of the least-squares fit is to minimize the square deviation, called χ^{2}, between the fitting function F (x) and the experimental data,

χ^{2} = ∑_{i=1}^{n} (\frac{F (x_{i}) − y_{i}}{Δy_{i}} )^{2} .  (6.6)

Individual deviations from experimental points are weighted with their inverse errors in order to promote contributions from the more precise measurements.

Minimization of χ^{2} with respect to the coefficiendt c_{k} in (6.5) is apparently equivalent to the least-squares problem (6.2) where

A_{ik} = \frac{f_{k}(x_{i})}{Δy_{i}} , b_{i} = \frac{y_{i}}{Δy_{i}} .  (6.7)

If QR = A is the QR-decomposition of the matrix A, the formal leastsquares solution to the fitting problem is

c = R^{−1}Q^{T} b .  (6.8)

However, in practice one has to back-substitute the system Rc = Q^{T} b.

6.4.1 Variances and correlations of fitting parameters

Suppose δy_{i} is a small deviation of the measured value of the physical observable at hand from its exact value. The corresponding deviation δc_{k} of the fitting coefficient is then given as

δc_{k} = ∑_{i} \frac{∂c_{k}}{∂y_{i}} δy_{i} .  (6.9)

In a good experiment the deviations δy_{i} are statistically independent and distributed normally with the standard deviations Δy_{i}. The deviations (6.9) are then also distributed normally with variances

⟨δc_{k}δc_{k}⟩ = ∑_{i} (\frac{∂c_{k}}{∂y_{i}} Δy_{i})^{2} = ∑_{i} (\frac{∂c_{k}}{∂b_{i}} )^{2} .  (6.10)

The standard errors in the fitting coefficients are then given as the square roots of variances,

Δc_{k} = \sqrt{⟨δc_{k}δc_{k}⟩} = \sqrt{∑_{i} (\frac{∂c_{k}}{∂b_{i}} )^{2}} .  (6.11)

The variances are diagonal elements of the covariance matrix, Σ, made of covariances,

Σ_{kq} ≡ ⟨δc_{k}δc_{q}⟩ = ∑_{i} \frac{∂c_{k}}{∂b_{i}} \frac{∂c_{q}}{∂b_{i}} .  (6.12)

Covariances ⟨δc_{k}δc_{q}⟩ are measures of to what extent the coefficients c_{k} and c_{q} change together if the measured values y_{i} are varied. The normalized covariances,

\frac{⟨δc_{k}δc_{q}⟩}{\sqrt{⟨δc_{k}δc_{k}⟩⟨δc_{q}δc_{q}⟩}}  (6.13)

are called correlations.

Using (6.12) and (6.8) the covariance matrix can be calculated as

Σ = (\frac{∂c}{∂b}) (\frac{∂c}{∂b})^{T} = R^{−1}(R^{−1})^{T} = (R^{T} R)^{−1} = (A^{T} A)^{−1} .  (6.14)

The square roots of the diagonal elements of this matrix provide the estimates of the errors Δc of the fitting coefficients,

Δc_{k} = \sqrt{Σ_{kk}} |_{k=1...m} ,  (6.15)

and the (normalized) off-diagonal elements provide the estimates of their correlations.

A C implementation of the ordinary least squares fit with QR-decomposition is shown in Table (6.1).

An illustration of a fit is shown on Figure (6.1) where a polynomial is fitted to a set of data.

60 CHAPTER 6. ORDINARY LEAST SQUARES

Table 6.1: C implementation of ordinary least squares fit with QRdecomposition.

#include <matrix . h>
#include <gc . h> // garbage c o l l e c t o r
#include <qr . h> // qr−decomposition
typedef struct { v e c t o r ∗ c ; matrix ∗S ; } l s f i t ;
l s f i t ∗ l s f i t a l l o c ( v e c t o r ∗x , v e c t o r ∗y , v e c t o r ∗dy ,
int nf , double f ( int k , double x ) )
{ int n=x−>s i z e , m=nf ;
matrix ∗A = m a t r i x a l l o c (n ,m) , ∗R = m a t r i x a l l o c (m,m) ;
v e c t o r ∗b = vector copy ( y ) ; v e c t o r d i v (b , dy ) ;
for ( int i =0; i<n ; i++)
{ double x i=v e c t o r g e t ( x , i ) , dyi=v e c t o r g e t ( dy , i ) ;
for ( int k=0;k<m; k++) m a t r i x s e t (A, i , k , f ( k , x i )/ dyi ) ;

}

qrdec (A,R) ;
l s f i t ∗ f i t = ( l s f i t ∗)GC MALLOC( sizeof ( l s f i t ) ) ;
f i t −>c = qrback (A,R, b ) ;
f i t −>S = q r i n v e r s e ( matrixT times matrix (R,R) ) ;
return f i t ;

}

Figure 6.1: Ordinary least squares fit of F_{c}(x) = c_{1} + c_{2}x + c_{3}x^{2} to a set of data. Shown are fits with optimal coefficiens c as well as with c + Δc and c − Δc.

Nonlinear equations

7.1 Introduction

Non-linear equations (or root-finding) is a problem of finding a set of n variables x = {x_{1}, …, x_{n}} which satisfy a system of n (generally) non-linear equations

f_{i}(x_{1}, ..., x_{n}) = 0 |_{i=1,...,n} .  (7.1)

In matrix notation the system is written as

f(x) = 0 ,  (7.2)

where f ≐ {f_{1}, …, f_{n}}.

In one-dimensional cases it is generally possible to plot the function in the region of interest and see whether the graph crosses the x-axis. One can then be sure the root exists and even figure out its approximate position to start one’s root-finding algorithm with. In multi-dimensions one generally does not know if the root exists at all, until it is found.

The root-finding algorithms generally proceed by iteration, starting from some approximate solution and making consecutive steps (hopefully) in the direction of the suspected root until some convergence criterion is satisfied. The procedure is generally not even guaranteed to converge unless starting from a point close enough to the sought root.

We shall only consider the multi-dimensional case here since i) the multidimensional root-finding is more difficult, and ii) the multi-dimensional routines can also be used in a one-dimensional case.

7.2 Newton’s method

Newton’s method (also reffered to as Newton-Raphson method, after Isaac Newton and Joseph Raphson) is a root-finding algorithm that uses the first term of the Taylor series of the functions f_{i} to linearise the system (7.1) in the vicinity of a suspected root. It is one of the oldest and best known methods and is a basis of a number of more refined methods.

Suppose that the point x = {x_{1}, …, x_{n}} is close to the root. The Newton’s algorithm tries to find the step Δx which would move the point towards the root, such that

f_{i}(x + Δx) = 0 |_{i=1,...,n} .  (7.3)

The first order Taylor expansion of (7.3) gives a system of linear equations,

f_{i}(x) + ∑_{k=1}^{n} \frac{∂f_{i}}{∂x_{k}} Δx_{k} = 0 |_{i=1,...,n} ,  (7.4)

or, in the matrix form,

JΔx = −f(x),  (7.5)

where f(x) ≐ {f_{0}(x), …, f_{n}(x)} and J is the matrix of partial derivatives^{1},

J_{ik} ≐ \frac{∂f_{i}}{∂x_{k}} ,  (7.6)

called the Jacobian matrix.

The solution Δx to the linear system (7.5)—called the Newton’s step— gives the approximate direction and the step-size towards the solution.

The Newton’s method converges quadratically if sufficiently close to the solution. Otherwise the full Newton’s step Δx might actually diverge from the solution. Therefore in practice a more conservative step λΔx with λ < 1 is usually taken. The strategy of finding the optimal λ is referred to as line search.

It is typically not worth the effort to find λ which minimizes ∥f(x + λΔx)∥ exactly, since Δx is only an approximate direction towards the root. Instead an inexact but quick minimization strategy is usually used, like the backtracking line search where one first attempts the full step, λ = 1, and then backtracks, λ ← λ/2, until the condition

∥f(x + λΔx)∥ < (1 − \frac{λ}{2}) ∥f(x)∥  (7.7)

is satisfied. If the condition is not satisfied for sufficiently small λ_{min} the step is taken with λ_{min} simply to step away from this diffictul place and try again. A typical algrorithm of the Newton’s method with backtracking line search and condition (7.7) is shown in Table 7.1.

\frac{∂f_{i}}{∂x_{k}} ≈ \frac{f_{i}(x_{1}, …, x_{k−1}, x_{k} + δx, x_{k+1}, …, x_{n}) − f_{i}(x_{1}, …, x_{k}, …, x_{n})}{δx}

1in practice if derivatives are not available analytically one uses finite differences with δx ≪ s where s is the typical scale of the problem at hand.

Table 7.1: Newton’s root-finding algorithm with backtracking.

repeat

solve JΔx = −f(x) for Δx

λ ← 1

while ( ∥f(x + λΔx)∥ > (1 − \frac{λ}{2} ) ∥f(x)∥ and λ > \frac{1}{64} ) do λ ← λ/2 x ← x + λΔx

until converged (e.g. ∥f(x)∥ < tolerance)

A somewhat more refined backtracking linesearch is based on an approximate minimization of the function

g(λ) ≐ \frac{1}{2}∥f(x + λΔx)∥^{2}  (7.8)

using interpolation. The values g(0) = \frac{1}{2}∥f(x)∥^{2} and g′(0) = −∥f(x)∥^{2} are already known (check this). If the previous step with certain λ_{trial} was rejected, we also have g(λ_{trial}). These three quantities allow to build a quadratic approximation,

g(λ) ≈ g(0) + g′(0)λ + cλ^{2} ,  (7.9)

where

c = \frac{g(λ_{trial}) − g(0) − g′(0)λ_{trial}}{λ_{trial}^{2}} .  (7.10)

The minimum of this approximation (determined by the condition g′(λ) = 0),

λ_{next} = −\frac{g′(0)}{2c} ,  (7.11)

becomes the next trial step-size.

The procedure is repeated recursively until either condition (7.7) is satisfied or the step becomes too small (in which case it is taken unconditionally in order to simply get away from the difficult place).

7.3 Broyden’s quasi-Newton method

The Newton’s method requires calculation of the Jacobian matrix at every iteration. This is generally an expensive operation. Quasi-Newton methods avoid calculation of the Jacobian matrix at the new point x + Δx, instead trying to use certain approximations, typically rank-1 updates.

Broyden’s algorithm estimates the Jacobian J + ΔJ at the point x + Δx using the finite-difference approximation,

(J + ΔJ)Δx = Δf ,  (7.12)

where Δf ≐ f(x + Δx) − f(x) and J is the Jacobian at the point x.

64 CHAPTER 7. NONLINEAR EQUATIONS

Table 7.2: C implementation of Newton’s algorithm with backtracking.

#include<matrix . h> // garbage c o l l e c t e d matrices and v e c t o r s
#include<qr . h>
v e c t o r ∗ newton
( v e c t o r ∗ f ( v e c t o r ∗x ) , v e c t o r ∗ startx , v e c t o r ∗ dx , double eps )

{

v e c t o r ∗x=vector copy ( s t a r t x ) , ∗ fx , ∗Dx ;
int n=x−>s i z e ;
matrix ∗A=m a t r i x a l l o c (n , n ) , ∗R=m a t r i x a l l o c (n , n ) ;

do {

fx=f ( x ) ; // : g e t s garbage c o l l e c t e d
for ( int j =0; j<n ; j++) // c a l c u l a t i n g the Jacobian matrix :

{

v e c t o r s e t ( x , j , v e c t o r g e t ( x , j )+ v e c t o r g e t ( dx , j ) ) ;
v e c t o r ∗ df=vector sum ( 1 , f ( x) , −1 , fx ) ;
for ( int i =0; i<n ; i++)
m a t r i x s e t (A, i , j , v e c t o r g e t ( df , i )/ v e c t o r g e t ( dx , j ) ) ;
v e c t o r s e t ( x , j , v e c t o r g e t ( x , j )− v e c t o r g e t ( dx , j ) ) ;

}

qrdec (A,R) ;
v e c t o r ∗ minusfx=vector copy ( fx ) ; v e c t o r s c a l e ( minusfx , −1);
Dx=qrback (A,R, minusfx ) ; // : Newton ’ s s t e p
double lambda=2, normfx=vector norm ( fx ) ;
v e c t o r ∗z , ∗ f z ;
do // back−t r a c k i n g l i n e −search :

{

lambda/=2; z=vector sum ( 1 , x , lambda , Dx ) ; f z=f ( z ) ;
}while ( vector norm ( f z )>(1−lambda /2)∗ normfx && lambda >0.02);
x=z ; fx=f z ;
}while ( vector norm (Dx)>vector norm ( dx ) && vector norm ( fx )>eps ) ;
return x ;

}

The matrix equation (7.12) is under-determined in more than one dimension as it contains only n equations to determine n^{2} matrix elements of ΔJ. Broyden suggested to choose ΔJ as a rank-1 update, linear in Δx,

ΔJ = c Δx^{T} ,  (7.13)

where the unknown vector c can be found by substituting (7.13) into (7.12), which gives

ΔJ = \frac{Δf − JΔx}{∥Δx∥^{2}}Δx^{T} .  (7.14)

In practice if one wanders too far from the point where J was first calculated the accuracy of the updates may decrease significantly. In such case one might need to recalculate J anew. For example, two successive steps with λ_{min} might be interpreted as a sign of accuracy loss in J and subsequently trigger its recalculation.

66 CHAPTER 7. NONLINEAR EQUATIONS

Minimization and optimization

8.1 Introduction

Minimization is a problem of finding the minimum (or the maximum) of a given, generally non-linear, real valued function f(x) (often called the objective function) of an n-dimensional argument x ≐ {x_{1}, …, x_{n}}.

Minimization is a simple case of a more general poblem — optimization — which includes finding best available values of some objective function within a given domain and subject to given constrains.

Minimization is related to root-finding as at the minimum all partial derivatives of the objective function vanish,

\frac{∂f}{∂x_{i}} = 0 .  (8.1)

One can alternatively solve this system of non-linear equations.

8.2 Local minimization

8.2.1 Newton’s methods

Newton’s method is based on the quadratic approxiamtion of the objective function f in the vicinity of the suspected minimum,

f(x + Δx) ≈ f(x) + ∇f(x)^{T} Δx + \frac{1}{2}Δx^{T} H(x)Δx ,  (8.2)

where the vector ∇f(x) is the gradient of the objective function at the point x,

∇f(x) ≐ [\frac{∂f(x)}{∂x_{i}}]_{i=1,...,n} ,  (8.3)

and H(x) is the Hessian matrix – a square matrix made of second-order partial derivatives of the objective function at the point x,

H(x) ≐ [\frac{∂^{2}f(x)}{∂x_{i}∂x_{j}}]_{i,j∈1,...,n} .  (8.4)

The minimum of this quadratic form as function of Δx is found at the point where its gradient with respect to Δx vanishes,

∇f(x) + H(x)Δx = 0 .  (8.5)

This gives an approximate step towards the minimum, called the Newton’s step,

Δx = −H(x)^{−1}∇f(x) .  (8.6)

The original Newton’s method is simply the iteration,

x_{k+1} = x_{k} − H(x_{k})^{−1}∇f(x_{k}) ,  (8.7)

where at each iteration the full Newton’s step is taken and the Hessian matrix is recalculated. In practice, instead of calculating H^{−1} one rather solves the linear equation (8.5).

Usually Newton’s method is modified to include a smaller step size, λΔx, with 0 < λ < 1. The step-size λ can be found by a backtracking algoritm similar to that in the Newton’s method for root-finding. One starts with λ = 1 and than backtracks, λ ← λ/2, until the Armijo condition,

f(x + λΔx) < f(x) + α∇f(x)λΔx ,  (8.8)

is satisfied. The parameter α can be chosen as small as 10^{−4}.

8.2.2 Quasi-Newton methods

Quasi-Newton methods are variations of the Newton’s method which attempt to avoid recalculation of the Hessian matrix at each iteration, trying instead certain updates based on the analysis of the gradient vectors. The update δH is usually chosen to satisfy the condition

∇f(x + s) = ∇f(x) + (H + δH)s ,  (8.9)

called secant equation, which is the Taylor expansion of the gradient. This secant condition is under-determined in more than one dimension as it consists of only n equations for the n^{2} unknown elements of the update δH. Various quasi-Newton methods use different choices for the form of the solution of the secant equation.

In quasi-Newton methods one often starts with a unity matrix as the zeroth approximation for the Hessian matrix and then applies the updates. In this case one can actually use the inverse Hessian matrix—thus avoiding the need to solve the linear equation at each iteration—and apply the updates directly to the inverse matrix.

Broydens’s update

The Broyden’s update is chosen in the form

δH = cs^{T} .  (8.10)

where the vector c is found from the condition (8.9). The resulting update is

H → H + \frac{y − Hs}{s^{T} s}s^{T} ,  (8.11)

where y ≐ ∇f(x + s) − ∇f(x). The update for the inverse matrix is given as^{1}

H^{−1} → H^{−1} + \frac{(s − H^{−1}y)s^{T} H^{−1}}{y^{T} H^{−1}s} .  (8.12)

SR1 update

The symmetric-rank-1 update (SR1) in chosen in the form

δH = uu^{T} ,  (8.13)

where the vector u is found from the condition (8.9). The resulting update is

H → H + \frac{(y − Hs)(y − Hs)^{T}}{(y − Hs)^{T} s} .  (8.14)

And for the inverse matrix

H^{−1} → H^{−1} + \frac{(s − H^{−1}y)(s − H^{−1}y)^{T}}{(s − H^{−1}y)^{T} y} .  (8.15)

Davidon–Fletcher–Powell (DFP) update

δH^{−1} = \frac{yy^{T}}{y^{T} s} − \frac{H^{−1}ss^{T} H^{−1}}{s^{T} H^{−1}s}  (8.16)

8.2.3 Downhill simplex method

The downhill simplex method (also called Nelder-Mead method or amoeba method) is a commonnly used minimization algorithm, where the minimum of a function in an n-dimensional space is found by transforming a simplex— a polytope with n+1 vertexes—according to the function values at the vertexes, moving it downhill until it converges towards the minimum.

(A + uv^{T})^{−1} = A^{−1} − \frac{A^{−1}uv^{T} A^{−1}}{1 + v^{T} A^{−1}u} .

1using the Sherman-Morrison formula,

The advantages of the downhill simplex method is its stability and the lack of use of derivatives. However, the convergence is realtively slow as compared to Newton’s methods.

In order to introduce the algorithm we need the following definitions:

• Simplex: a figure (polytope) represented by n+1 points, called vertexes, {p_{1}, …, p_{n+1}} (where each point p_{k} is an n-dimensional vector).

• Highest point: the vertex, p_{hi}, with the highest value of the function: f(p_{hi}) = max_{k} f(p_{k}).

• Lowest point: the vertex, p_{lo}, with the lowest value of the function: f(p_{lo}) = min_{k} f(p_{k}).

• Centroid: the center of gravity of all points, except for the highest: p_{ce} = \frac{1}{n} ∑_{(k≠hi)} p_{k}

The simplex is moved downhill by a combination of the following elementary operations:

1. Reflection: the highest point is reflected against the centroid, p_{hi} → p_{re} = p_{ce} + (p_{ce} − p_{hi}).

2. Expansion: the highest point reflects and then doubles its distance from the centroid, p_{hi} → p_{ex} = p_{ce} + 2(p_{ce} − p_{hi}).

3. Contraction: the highest point halves its distance from the centroid, p_{hi} → p_{co} = p_{ce} + \frac{1}{2}(p_{hi} − p_{ce}).

4. Reduction: all points, except for the lowest, move towards the lowest points halving the distance. p_{k≠lo} → \frac{1}{2}(p_{k} + p_{lo}).

A C-implementation of the simplex operations is listed in Table 8.4.

Finally, Table 8.1 shows one possible algorithm for the downhill simplex method.

8.3 Global optimization

Global optimization is the problem of locating (a good approximation to) the global minimum of a given objective function in a search space large enough to prohibit exhaustive enumeration.

When only a small sub-space of the search space can be realistically sampled the stochastic methods usually come to the fore.

A good local minimizer converges to the nearest local miminum relatively fast, so one possible global minimizer can be constructed by simply starting the local miminizer from different random starting points.

In the following several popular global minimization algorithms are shortly described.

Table 8.1: Downhill simplex (Nelder-Mead) algorithm

REPEAT :
f i n d highest , lowest , and c e n t r o i d p o i n t s o f the simplex
try r e f l e c t i o n
IF f(reflected) < f(lowest) :
try expansion
IF f(expanded) < f(reflected) :
accept expansion
ELSE :
accept r e f l e c t i o n
ELSE :
IF f(reflected) < f(highest) :
accept r e f l e c t i o n
ELSE :
try c o n t r a c t i o n
IF f(contracted) < f(highest) :
accept c o n t r a c t i o n
ELSE :
do r e d u c t i o n
UNTIL converged ( e . g . s i z e ( simplex)< t o l e r a n c e )

8.3.1 Simulated annealing

Simulated annealing is a stochastic metaheuristic algorithm for global minimization. The name and inspiration come from annealing — heating up and cooling slowly — in material science. The slow cooling allows a piece of material to reach a state with ”lowest energy”.

The objective function in the space of states is interpreted as some sort of potential energy and the states — the points in the search space — are interpeted as physical states of some physical system. The system attempts to make transitions from its current state to some randomly sampled nearest states with the goal to eventually reach the state with minimal energy – the global minimum.

The system is attached to a thermal resevoir with certain temperature. Each time the energy of the system is measured the reservoir supplies it with a random amount of thermal energy sampled from the Boltzmann distribution,

P (E) = T e^{−E/T} .  (8.17)

If the temperature equals zero the system can only make transitions to the neighboring states with lower potential energy. In this case the algorithm turns merely into a local minimizer with random sampling.

If temperature is finite the system is able to climb up the ridges of the potential energy — about as high as the current temperature — and thus escape from local minima and hopefully eventually reach the global minimum.

Table 8.2: Simulated annealing algorithm

s t a t e ← s t a r t s t a t e
T ← s t a r t t e m p e r a t u r e
energy ← E( s t a t e ) + random thermal energy (T)
REPEAT :
new state ← neighbour ( s t a t e )
new energy ← E( new state ) + random thermal energy (T)
IF new energy < energy :
s t a t e ← new state
energy ← new energy
r e d u c e t e m p e r a t u r e a c c o r d i n g t o s c h e d u l e (T)
UNTIL terminated

One typically starts the simulation with some finite temperature on the order of the height of the typical hills on the potential energy surface, letting the system to wander almost unhindered around the landscape with a good chance to locale if not the best then at least a good enough minimum. The temperature is then slowly reduced following some annealing schedule which may be supplied by the user but must end with T = 0 towards the end of the alloted time budget.

Table 8.2 lists one possible variant of the algorithm.

The function neigbour(state) should return a randomly chosen neighbour of the given state.

Downhill simplex method can incorporate simulated annealing by adding the stochastic thermal energy to the values of the objective function at the verices.

8.3.2 Quantum annealing

Quantum annealing is a general global minimization algorithm which — like simulated annealing — also allows the search path to escape from local minima. However instead of the thermal jumps over the potential barriers quantum annealing allows the system to tunnel through the barriers.

In its simple incarnation the quantum annealing algorithm allows the system to attempt transitions not only to the nearest states but also to distant states within certain ”tunneling distance” from the current state. The transition is accepted only if it reduces the potential energy of the system.

At the beginning of the minimization prcedure the tunnelling distance is large — on the order of the size of the region where the global minimum is suspected to be localed — allowing the system to explore the region. The tunneling distance is then slowly reduced according to a schedule such that by the end of the alloted time the tunnelling distance reduces to zero at which point the system hopefully is in the state with minimal energy.

Table 8.3: Quantum annealing algorithm

s t a t e ← s t a r t s t a t e
energy ← E( s t a t e )
R ← s t a r t r a d i u s
REPEAT :
new state ← n e i g h b o u r w i t h i n r a d i u s ( state ,R)
new energy ← E( new state )
IF new energy < energy :
s t a t e ← new state
energy ← new energy
r e d u c e r a d i u s a c c o r d i n g t o s c h e d u l e (R)
UNTIL terminated

8.3.3 Evolutionary algorithms

Unlike annealing algorithms, which follow the motion of only one point in the search space, the evolutionary algorithms typically follow a set of points called a population of individuals. A bit like the downhill simplex method which follows the motion of a set of points – the simplex.

The population evolves towards more fit individuals where fitness is understood in the sence of minimizing the objective functions. The parameters of the individuals (for example, the coordinates of the points in multi-dimentional minimization of a continuous objective function) are called genes.

The algorithm proceeds iteratively with the population in each iteration called a generation. In each generation the fitness of each individual — typically, the value of the objective function — is evaluated and the new generation is generated stochastically from the gene pool of the current generation through crossovers and mutations such that the genes of more fit individuals have a better chance of propagating into the next generation.

Each new individual in the next generation is produced from a pair of ”parent” individuals of the current generation. The use of two ”parents” is biologically inspired, in practice more than two ”parents” can be used as well. The parents for a new individual are selected from the individuals of the current generation through a fitness based stochastic process where fitter individuals are more likely to be selected.

The ”child” individual shares many characteristics of its ”parents”. In the simplest case the ”child” may get its genes by simply averaging the genes of its parents. Then a certain amount of mutations — random changes in the genes — are added to the ”child’s” genes.

Generation of ”children” continues until the population of the new generation reaches the appropriate size after which the iteration repeats itself.

The algorithm is terminated when the fitness level of the population is deemed sufficient or when the allocated budget is exhausted.

74 CHAPTER 8. MINIMIZATION AND OPTIMIZATION

8.4 Implementation in C

8.4. IMPLEMENTATION IN C 75

Table 8.4: C implementation of simplex operations

#include <math . h>
#define float double
void r e f l e c t i o n
( float ∗ highest , float ∗ centroid , int dim , float ∗ r e f l e c t e d ){
for ( int i =0; i<dim ; i++) r e f l e c t e d [ i ]=2∗ c e n t r o i d [ i ]− h i g h e s t [ i ] ;

}

void expansion
( float ∗ highest , float ∗ centroid , int dim , float ∗ expanded ) {
for ( int i =0; i<dim ; i++) expanded [ i ]=3∗ c e n t r o i d [ i ]−2∗ h i g h e s t [ i ] ;

}

void c o n t r a c t i o n
( float ∗ highest , float ∗ centroid , int dim , float ∗ contracted ){
for ( int i =0; i<dim ; i++)
contracted [ i ]=0.5∗ c e n t r o i d [ i ]+0.5∗ h i g h e s t [ i ] ;

}

void r e d u c t i o n ( float ∗∗ simplex , int dim , int l o ){
for ( int k=0;k<dim+1;k++) i f ( k!= l o ) for ( int i =0; i<dim ; i++)
simplex [ k ] [ i ]=0.5∗( simplex [ k ] [ i ]+ simplex [ l o ] [ i ] ) ;

}

float d i s t a n c e ( float ∗ a , float ∗ b , int dim ){
float s =0; for ( int i =0; i<dim ; i++) s+=pow( b [ i ]−a [ i ] , 2 ) ;
return s q r t ( s ) ;

}

float s i z e ( float ∗∗ simplex , int dim ){
float s =0; for ( int k=1;k<dim+1;k++){
float d i s t=d i s t a n c e ( simplex [ 0 ] , simplex [ k ] , dim ) ;
i f ( d i s t >s ) s=d i s t ; }
return s ;

}

#define float double
void simplex update ( float ∗∗ simplex , float ∗ f v a l u e s , int d ,
int ∗ hi , int ∗ lo , float ∗ c e n t r o i d )

{

∗ hi =0; float h i g h e s t=f v a l u e s [ 0 ] ;
∗ l o =0; float lowest =f v a l u e s [ 0 ] ;
for ( int k=1;k<d+1;k++) {
float next=f v a l u e s [ k ] ;
i f ( next>h i g h e s t ){ h i g h e s t=next ; ∗ hi=k ; }
i f ( next<lowest ) { lowest=next ; ∗ l o=k ; }
}
for ( int i =0; i<d ; i++) {
float s =0; for ( int k=0;k<d+1;k++) i f ( k!=∗ hi ) s+=simplex [ k ] [ i ] ;
c e n t r o i d [ i ]= s /d ; }

}

void s i m p l e x i n i t i a t e (
float fun ( float ∗) , float ∗∗ simplex , float ∗ f v a l u e s , int d ,
int ∗ hi , int ∗ lo , float ∗ c e n t r o i d )

{

for ( int k=0;k<d+1;k++) f v a l u e s [ k]= fun ( simplex [ k ] ) ;
simplex update ( simplex , f v a l u e s , d , hi , lo , c e n t r o i d ) ;

}

76 CHAPTER 8. MINIMIZATION AND OPTIMIZATION

Table 8.5: C implementation of downhill simplex algorithm

#include <simplex . h>
#include <s t d l i b . h>
#include <s t d i o . h>
int d o w n h i l l s i m p l e x (
double F( double ∗) , double∗∗ simplex , int d , double s i m p l e x s i z e g o a l )

{

int hi , lo , k=0; double c e n t r o i d [ d ] , F value [ d+1] , p1 [ d ] , p2 [ d ] ;
s i m p l e x i n i t i a t e (F, simplex , F value , d,& hi ,& lo , c e n t r o i d ) ;
while ( s i z e ( simplex , d)> s i m p l e x s i z e g o a l ){
simplex update ( simplex , F value , d,& hi ,& lo , c e n t r o i d ) ;
r e f l e c t i o n ( simplex [ hi ] , centroid , d , p1 ) ; double f r e=F( p1 ) ;
i f ( f r e <F value [ l o ] ) { // r e f l e c t i o n l o o k s good : t r y expansion
expansion ( simplex [ hi ] , centroid , d , p2 ) ; double f e x=F( p2 ) ;
i f ( f e x <f r e ){ // accept expansion
for ( int i =0; i<d ; i++)simplex [ hi ] [ i ]=p2 [ i ] ; F value [ hi ]= f e x ; }
else { // r e j e c t expansion and accept r e f l e c t i o n
for ( int i =0; i<d ; i++)simplex [ hi ] [ i ]=p1 [ i ] ; F value [ hi ]= f r e ; } }
else { // r e f l e c t i o n wasn ’ t good
i f ( f r e <F value [ hi ] ) { // ok , accept r e f l e c t i o n
for ( int i =0; i<d ; i++)simplex [ hi ] [ i ]=p1 [ i ] ; F value [ hi ]= f r e ; }
else { // t r y c o n t r a c t i o n
c o n t r a c t i o n ( simplex [ hi ] , centroid , d , p1 ) ; double f c o=F( p1 ) ;
i f ( f c o <F value [ hi ] ) { // accept c o n t r a c t i o n
for ( int i =0; i<d ; i++)simplex [ hi ] [ i ]=p1 [ i ] ; F value [ hi ]= f c o ; }
else { // do r e d u c t i o n
r e d u c t i o n ( simplex , d , l o ) ;
s i m p l e x i n i t i a t e (F, simplex , F value , d,& hi ,& lo , c e n t r o i d ) ; } } }
k++;} return k ;

}

Eigenvalues and eigenvectors

9.1 Introduction

A non-zero column-vector v is called an eigenvector of a matrix A with an eigenvalue λ, if

Av = λv .  (9.1)

If an n × n matrix A is real and symmetric, A^{T} = A, then it has n real eigenvalues λ_{1}, …, λ_{n}, and its (orthogonalized) eigenvectors V = {v_{1}, …, v_{n}} form a full basis,

V V^{T} = V^{T} V = 1 ,  (9.2)

in which the matrix is diagonal,

V^{T} AV = [\begin{matrix}λ_{1} & 0 & ⋯ & 0 \\ 0 & λ_{2} & & ⋮ \\ ⋮ & & ... \\ 0 & ⋯ & & λ_{n} \end{matrix}] .  (9.3)

Matrix diagonalization means finding all eigenvalues and (optionally) eigenvectors of a matrix.

Eigenvalues and eigenvectors enjoy a multitude of applications in different branches of science and technology.

9.2 Similarity transformations

Orthogonal transformations,

A → Q^{T} AQ ,  (9.4)

where Q^{T} Q = 1, and, generally, similarity transformations,

A → S^{−1}AS ,  (9.5)

preserve eigenvalues and eigenvectors. Therefore one of the strategies to diagonalize a matrix is to apply a sequence of similarity transformations (also called rotations) which (iteratively) turn the matrix into diagonal form.

9.2.1 Jacobi eigenvalue algorithm

Jacobi eigenvalue algorithm is an iterative method to calculate the eigenvalues and eigenvectors of a real symmetric matrix by a sequence of Jacobi rotations.

Jacobi rotation is an orthogonal transformation which zeroes a pair of the off-diagonal elements of a (real symmetric) matrix A,

A → A′ = J(p, q)^{T} AJ(p, q) : A′_{pq} = A′_{qp} = 0 .  (9.6)

The orthogonal matrix J(p, q) which eliminates the element A_{pq} is called the Jacobi rotation matrix. It is equal identity matrix except for the four elements with indices pp, pq, qp, and qq,

J(p, q) = [\begin{matrix}1 \\ & ... & & & & 0 \\ & & cos ϕ & ⋯ & sin ϕ \\ & & ⋮ & ... & ⋮ \\ & & − sin ϕ & ⋯ & cos ϕ \\ & 0 & & & & ... \\ & & & & & & 1 \end{matrix}] \begin{matrix}← row p \\ ← row q \end{matrix} .  (9.7)

Or explicitly,

J(p, q)_{ij} = δ_{ij} ∀ ij ∉ {pq, qp, pp, qq} ;
J(p, q)_{pp} = cos ϕ = J(p, q)_{qq} ;
J(p, q)_{pq} = sin ϕ = −J(p, q)_{qp} .  (9.8)

After a Jacobi rotation, A → A′ = J^{T} AJ, the matrix elements of A′ become

A′_{ij} = A_{ij} ∀ i ≠ p, q ∧ j ≠ p, q
A′_{pi} = A′_{ip} = cA_{pi} − sA_{qi} ∀ i ≠ p, q ;
A′_{qi} = A′_{iq} = sA_{pi} + cA_{qi} ∀ i ≠ p, q ;
A′_{pp} = c^{2}A_{pp} − 2scA_{pq} + s^{2}A_{qq} ;
A′_{qq} = s^{2}A_{pp} + 2scA_{pq} + c^{2}A_{qq} ;
A′_{pq} = A′_{qp} = sc(A_{pp} − A_{qq}) + (c^{2} − s^{2})A_{pq} ,  (9.9)

where c ≡ cos ϕ, s ≡ sin ϕ. The angle ϕ is chosen such that after rotation the matrix element A′_{pq} is zeroed,

cot(2ϕ) = \frac{A_{qq} − A_{pp}}{2A_{pq}} ⇒ A′_{pq} = 0 .  (9.10)

A side effect of zeroing a given off-diagonal element A_{pq} by a Jacobi rotation is that other off-diagonal elements are changed. Namely, the elements of the rows and columns with indices equal to p and q. However, after the Jacobi rotation the sum of squares of all off-diagonal elements is reduced. The algorithm repeatedly performs rotations until the off-diagonal elements become sufficiently small.

The convergence of the Jacobi method can be proved for two strategies for choosing the order in which the elements are zeroed:

1. Classical method: with each rotation the largest of the remaining offdiagonal elements is zeroed.

2. Cyclic method: the off-diagonal elements are zeroed in strict order, e.g. row after row.

Although the classical method allows the least number of rotations, it is typically slower than the cyclic method since searching for the largest element is an O(n^{2}) operation. The count can be reduced by keeping an additional array with indexes of the largest elements in each row. Updating this array after each rotation is only an O(n) operation.

A sweep is a sequence of Jacobi rotations applied to all non-diagonal elements. Typically the method converges after a small number of sweeps. The operation count is O(n) for a Jacobi rotation and O(n^{3}) for a sweep.

The typical convergence criterion is that the sum of absolute values of the off-diagonal elements is small, ∑_{i<j} |A_{ij}| < ϵ, where ϵ is the required accuracy. Other criteria can also be used, like the largest off-diagonal element is small, max |A_{i<j}| < ϵ, or the diagonal elements have not changed after a sweep.

The eigenvectors can be calculated as V = 1J_{1}J_{2}..., where J_{i} are the successive Jacobi matrices. At each stage the transformation is

V_{ij} → V_{ij} , j ≠ p, q
V_{ip} → cV_{ip} − sV_{iq}  (9.11)
V_{iq} → sV_{ip} + cV_{iq}

Alternatively, if only one (or few) eigenvector v_{k} is needed, one can instead solve the (singular) system (A − λ_{k})v = 0.

QR/QL algorithm

An orthogonal transformation of a real symmetric matrix, A → Q^{T} AQ = RQ, where Q is from the QR-decomposition of A, partly turns the matrix A into diagonal form. Successive iterations eventually make it diagonal. If there are degenerate eigenvalues there will be a corresponding block-diagonal sub-matrix. For convergence properties it is of advantage to use shifts: instead of QR[A] we do QR[A − s1] and then A → RQ + s1. The shift s can be chosen as A_{nn}. As soon as an eigenvalue is found the matrix is deflated, that is the corresponding row and column are crossed out.

Accumulating the successive transformation matrices Q_{i} into the total matrix Q = Q_{1} … Q_{N}, such that Q^{T} AQ = Λ, gives the eigenvectors as columns of the Q matrix.

If only one (or few) eigenvector v_{k} is needed one can instead solve the (singular) system (A − λ_{k})v = 0.

Tridiagonalization.

Each iteration of the QR/QL algorithm is an O(n^{3}) operation. On a tridiagonal matrix it is only O(n). Therefore the effective strategy is first to make the matrix tridiagonal and then apply the QR/QL algorithm. Tridiagonalization of a matrix is a non-iterative operation with a fixed number of steps.

9.3 Eigenvalues of updated matrix

In practice it happens quite often that the (symmetric) matrix A to be diagonalized is given in the form of a diagonal matrix, D, plus an update matrix, W,

A = D + W ,  (9.12)

where the update W is a simpler — in a certain sense — matrix which allows a more efficient calculation of the updated eigenvalues, as compared to general diagonalization algorithms.

The most common updates are

• symmetric rank-1 update,

W = uu^{T} ,  (9.13)

where u is a columnt-vector;

• symmetric rank-2 update,

W = uv^{T} + vu^{T} ;  (9.14)

• symmetric row/column update – a special case of rank-2 update,

W = [\begin{matrix}0 & … & u_{1} & … & 0 \\ ⋮ & ... & ⋮ & ... & ⋮ \\ u_{1} & … & u_{p} & … & u_{n} \\ ⋮ & ... & ⋮ & ... & ⋮ \\ 0 & … & u_{n} & … & 0 \end{matrix}] ≡ e(p)u^{T} + ue(p)^{T} ,  (9.15)

where e(p) is the unit vector in the p-direction.

9.3.1 Rank-1 update

We assume that the size-n real symmetric matrix A to be diagonalized is given in the form of a diagonal matrix plus a rank-1 update,

A = D + σuu^{T} ,  (9.16)

where D is a diagonal matrix with diagonal elements {d_{1}, …, d_{n}} and u is a given vector. The diagonalization of such matrix can be done in O(m^{2}) operations, where m ≤ n is the number of non-zero elements in the update vector u, as compared to O(n^{3}) operations for a general diagonalization [7].

The eigenvalue equation for the updated matrix reads

(D + σuu^{T} ) q = λq ,  (9.17)

where λ is an eigenvalue and q is the corresponding eigenvector. The equation can be rewritten as

(D − λ1) q + σuu^{T} q = 0 .  (9.18)

Multiplying from the left with u^{T} (D − λ1)^{−1} gives

u^{T} q + u^{T} (D − λ1)^{−1} σuu^{T} q = 0 .  (9.19)

Finally, dividing by u^{T} q leads to the (scalar) secular equation (or characteristic equation) in λ,

1 + ∑_{i=1}^{m} \frac{σu_{i}^{2}}{d_{i} − λ} = 0 ,  (9.20)

where the summation index counts the m non-zero components of the update vector u. The m roots of this equation determine the (updated) eigenvalues^{1}.

Finding a root of a rational function requires an iterative technique, such as the Newton-Raphson method. Therefore diagonalization of an updated matrix is still an iterative procedure. However, each root can be found in O(1) iterations, each iteration requiring O(m) operations. Therefore the iterative part of this algorithm needs O(m^{2}) operations.

Finding roots of this particular secular equation can be simplified by utilizing the fact that its roots are bounded by the eigenvalues d_{i} of the matrix D. Indeed if we denote the roots as λ_{1}, λ_{2}, …, λ_{n} and assume that λ_{i} ≤ λ_{i+1} and d_{i} ≤ d_{i+1}, it can be shown that

1. if σ ≥ 0,

d_{i} ≤ λ_{i} ≤ d_{i+1} ,  i = 1, …, n − 1 ,  (9.21)
d_{n} ≤ λ_{n} ≤ d_{n} + σu^{T} u ;  (9.22)

1Multiplying this equation by ∏_{i=1}^{m}(d_{i} − λ) leads to an equivalent polynomial equation of the order m, which has exactly m roots.

2. if σ ≤ 0,

d_{i−1} ≤ λ_{i} ≤ d_{i} ,  i = 2, …, n ,  (9.23)
d_{1} + σu^{T} u ≤ λ_{1} ≤ d_{1} .  (9.24)

9.3.2 Symmetric row/column update

The matrix A to be diagonalized is given in the form

A = D + e(p)u^{T} + ue(p)^{T} = [\begin{matrix}d_{1} & … & u_{1} & … & 0 \\ ⋮ & ... & ⋮ & ... & ⋮ \\ u_{1} & … & d_{p} & … & u_{n} \\ ⋮ & ... & ⋮ & ... & ⋮ \\ 0 & … & u_{n} & … & d_{n} \end{matrix}] ,  (9.25)

where D is a diagonal matrix with diagonal elements {d_{i}|i = 1, …, n}, e(p) is the unit vector in the p-direction, and u is a given update vector where the p-th element can be assumed to equal zero, u_{p} = 0, without loss of generality. Indeed, if the element is not zero, one can simply redefine d_{p} → d_{p} + 2u_{p}, u_{p} → 0.

The eigenvalue equation for matrix A is given as

(D − λ)x + e(p)u^{T} x + ue(p)^{T} x = 0 ,  (9.26)

where x is an eigenvector and λ is the corresponding eigenvalue. The component number p of this vector-equation reads

(d_{p} − λ)x_{p} + u^{T} x = 0 ,  (9.27)

while the component number k ≠ p reads

(d_{k} − λ)x_{k} + u_{k}x_{p} = 0 ,  (9.28)

Dividing the last equation by (d_{k} − λ), multiplying from the left with ∑_{k=1}^{n} u_{k}, substituting u^{T} x using equation (9.27) and dividing by x_{p} gives the secular equation,

−(d_{p} − λ) + ∑_{k≠p}^{n} \frac{u_{k}^{2}}{d_{k} − λ} = 0 ,  (9.29)

which determines the updated eigenvalues.

(9.30)

9.3.3 Symmetric rank-2 update

A symmetric rank-2 update can be represented as two consecutive rank-1 updates,

uv^{T} + vu^{T} = aa^{T} − bb^{T} ,  (9.31)

where

a = \frac{1}{\sqrt{2}}(u + v) , b = \frac{1}{\sqrt{2}}(u − v) .  (9.32)

The eigenvalues can then be found by applying the rank-1 update method twice.

84 CHAPTER 9. EIGENVALUES AND EIGENVECTORS

Table 9.1: Jacobi diagonalization in C. The header matrix.h is supposed to define the self-explanatory functions vector get, vector set, matrix get,

matrix set, and matrix set identity.

#include<math . h>
#include<matrix . h>
int j a c o b i ( matrix ∗ A, v e c t o r ∗ e , matrix ∗ V){
// Jacobi d i a g o n a l i z a t i o n .
// Upper t r i a n g l e o f A i s d e s t r o y e d .
// e and V accumulate e i g e n v a l u e s and e i g e n v e c t o r s
int changed , r o t a t i o n s =0, n=A−>s i z e 1 ;
for ( int i =0; i<n ; i++)v e c t o r s e t ( e , i , matrix get (A, i , i ) ) ;
m a t r i x s e t i d e n t i t y (V) ;
do{
changed =0; int p , q ;
for ( p=0;p<n ; p++)for ( q=p+1;q<n ; q++){
double app=v e c t o r g e t ( e , p ) ;
double aqq=v e c t o r g e t ( e , q ) ;
double apq=matrix get (A, p , q ) ;
double phi =0.5∗ atan2 (2∗ apq , aqq−app ) ;
double c = cos ( phi ) ;
double s = s i n ( phi ) ;
double app1=c ∗ c ∗app−2∗s ∗ c ∗apq+s ∗ s ∗aqq ;
double aqq1=s ∗ s ∗app+2∗s ∗ c ∗apq+c ∗ c ∗aqq ;
i f ( app1!=app | | aqq1!=aqq ){
changed =1; r o t a t i o n s ++;
v e c t o r s e t ( e , p , app1 ) ;
v e c t o r s e t ( e , q , aqq1 ) ;
m a t r i x s e t (A, p , q , 0 . 0 ) ;
for ( int i =0; i<p ; i ++){
double aip=matrix get (A, i , p ) ;
double aiq=matrix get (A, i , q ) ;
m a t r i x s e t (A, i , p , c ∗ aip−s ∗ aiq ) ;
m a t r i x s e t (A, i , q , c ∗ aiq+s ∗ aip ) ;

}

for ( int i=p+1; i<q ; i ++){
double api=matrix get (A, p , i ) ;
double aiq=matrix get (A, i , q ) ;
m a t r i x s e t (A, p , i , c ∗ api−s ∗ aiq ) ;
m a t r i x s e t (A, i , q , c ∗ aiq+s ∗ api ) ;

}

for ( int i=q+1; i<n ; i ++){
double api=matrix get (A, p , i ) ;
double aqi=matrix get (A, q , i ) ;
m a t r i x s e t (A, p , i , c ∗ api−s ∗ aqi ) ;
m a t r i x s e t (A, q , i , c ∗ aqi+s ∗ api ) ;

}

for ( int i =0; i<n ; i ++){
double vip=matrix get (V, i , p ) ;
double viq=matrix get (V, i , q ) ;
m a t r i x s e t (V, i , p , c ∗ vip−s ∗ viq ) ;
m a t r i x s e t (V, i , q , c ∗ viq+s ∗ vip ) ;

} }

}

}while ( changed ! = 0 ) ;
return r o t a t i o n s ;

}

Power iteration methods and Krylov subspaces

10.1 Introduction

Power method is an iterative method to calculate an eigenvalue and the corresponding eigenvector of a matrix A using the power iteration

x_{i+1} = Ax_{i} .  (10.1)

The iteration converges to the eigenvector with the largest eigenvalue. The eigenvalue can be estimated using the Rayleigh quotient

λ[x_{i}] = \frac{x_{i}^{T} Ax_{i}}{x_{i}^{T} x_{i}} = \frac{x_{i+1}^{T}x_{i}}{x_{i}^{T} x_{i}} .  (10.2)

The inverse power iteration with the inverse matrix,

x_{i+1} = A^{−1}x_{i} ,  (10.3)

converges to the smallest eigenvalue of matrix A.

Alternatively, the shifted inverse iteration,

x_{i+1} = (A − s1)^{−1}x_{i} ,  (10.4)

where 1 signifies a unity matrix of the same size as A, converges to the eigenvalue closest to the given number s.

The inverse iteration method is a refinement of the inverse power method where the trick is not to invert the matrix in (10.4) but rather solve the linear system

(A − s1)x_{i+1} = x_{i}  (10.5)

using e.g. QR-decomposition.

86CHAPTER 10. POWER ITERATION METHODS AND KRYLOV SUBSPACES

The better approximation s to the sought eigenvalue is chosen, the faster convergence one gets. However, incorrect choice of s can lead to slow convergence or to the convergence of a different eigenvector. Usually in practice the method is used when good approximation for the eigenvalue is known, and hence one needs only few (quite often just one) iteration.

One can update the estimate for the eigenvalue using the Rayleigh quotient λ[x_{i}] after each iteration and get faster convergence for the price of O(n^{3}) operations per QR-decomposition; or one can instead make more iterations (with O(n^{2}) operations per iteration) using the same matrix (A−s1). The optimal strategy is probably an update after several iterations.

10.2 Krylov subspaces

When calculating an eigenvalue of a matrix A using the power method, one starts with an initial random vector b and then computes iteratively the sequence Ab, A^{2}b, …, A^{n−1}b normalising and storing the result in b on each iteration. The sequence converges to the eigenvector of the largest eigenvalue of A.

The set of vectors

K_{n} = {b, Ab, A^{2}b, …, A^{n−1}b} ,  (10.6)

where n < rank(A), is called the order-n Krylov matrix, and the subspace spanned by these vectors is called the order-n Krylov subspace. The vectors are not orthogonal but can be made so e.g. by Gram-Schmidt orthogonalisation.

For the same reason that A^{n−1}b approximates the dominant eigenvector one can expect that the other orthogonalised vectors approximate the eigenvectors of the n largest eigenvalues.

Krylov subspaces are the basis of several successful iterative methods in numerical linear algebra, in particular: Arnoldi and Lanczos methods for finding one (or a few) eigenvalues of a matrix; and GMRES (Generalised Minimum RESidual) method for solving systems of linear equations.

These methods are particularly suitable for large sparse matrices as they avoid matrix-matrix operations but rather multiply vectors by matrices and work with the resulting vectors and matrices in Krylov subspaces of modest sizes.

10.3 Arnoldi iteration

Arnoldi iteration is an algorithm where the order-n orthogonalised Krylov matrix Q_{n} for a given matrix A is built using stabilised Gram-Schmidt process:

start with a set Q = {q^{1}} consisting of one random normalised vector q^{1};
repeat for k = 2 to n :
   make a new vector q^{k} = Aq^{k−1}
   orthogonalise q^{k} to all vectors q^{i} ∈ Q storing q^{i}^{†}q^{k} → h_{i,k−1}
   normalise q^{k} storing ∥q^{k}∥ → h_{k,k−1}
   add q^{k} to the set Q

By construction the matrix H_{n} made of the elements h_{jk} is an upper Hessenberg matrix,

H_{n} = [\begin{matrix}h_{1,1} & h_{1,2} & h_{1,3} & ⋯ & h_{1,n} \\ h_{2,1} & h_{2,2} & h_{2,3} & ⋯ & h_{2,n} \\ 0 & h_{3,2} & h_{3,3} & ⋯ & h_{3,n} \\ ⋮ & ... & ... & ... & ⋮ \\ 0 & ⋯ & 0 & h_{n,n−1} & h_{n,n} \end{matrix}] ,  (10.7)

which is a partial orthogonal reduction of A into Hessenberg form,

H_{n} = Q_{n}^{†}AQ_{n} .  (10.8)

The matrix H_{n} can be viewed as a representation of A in the Krylov subspace K_{n}. The eigenvalues and eigenvectors of the matrix H_{n} approximate the largest eigenvalues of matrix A.

Since H_{n} is a Hessenberg matrix of modest size its eigenvalues can be relatively easily computed with standard algorithms.

In practice if the size n of the Krylov subspace becomes too large the method is restarted.

10.4 Lanczos iteration

Lanczos iteration is Arnoldi iteration for Hermitian matrices, in which case the Hessenberg matrix H_{n} of Arnoldi method becomes a tridiagonal matrix T_{n}.

The Lanczos algorithm thus reduces the original hermitian N ×N matrix A into a smaller n × n tridiagonal matrix T_{n} by an orthogonal projection onto the order-n Krylov subspace. The eigenvalues and eigenvectors of a tridiagonal matrix of a modest size can be easily found by e.g. the QRdiagonalisation method.

In practice the Lanczos method is not very stable due to round-off errors leading to quick loss of orthogonality. The eigenvalues of the resulting tridiagonal matrix may then not be a good approximation to the original matrix. Library implementations fight the stability issues by trying to prevent the loss of orthogonality and/or to recover the orthogonality after the basis is generated.

10.5 Generalised minimum residual (GMRES)

GMRES is an iterative method for the numerical solution of a system of linear equations,

Ax = b ,  (10.9)

where the exact solution x is approximated by the vector x_{n} ∈ K_{n} that minimises the residual Ax_{n} − b in the Krylov subspace K_{n} of matrix A,

x ≈ x_{n} ← min_{x∈K_{n}} ∥Ax − b∥ .  (10.10)

The vector x_{n} ∈ K_{n} can be represented as x_{n} = Q_{n}y_{n} where Q_{n} is the projector on the space K_{n} and y_{n} is an n-dimensional vector. Substituting this into (10.9) gives

AQ_{n}y_{n} = b ,  (10.11)

which can be solved by least-squares method.

Fast Fourier transform

Fast Fourier transform (FFT) is an efficient algorithm to compute the discrete Fourier transform (DFT).

Computing DFT of N points in the naive way, using the definition, takes O(N^{2}) arithmetic operations, while an FFT can compute the same result in only O(N log N) operations. The difference in speed can be substantial, especially for large data sets. This improvement made many DFT-based algorithms practical.

Since the inverse of a DFT is also a DFT any FFT algorithm can be used in for the inverse DFT as well.

The most well known FFT algorithms, like the Cooley-Tukey algorithm, depend upon the factorization of N. However, there are FFTs with O(N log N) complexity for all N, even for prime N.

11.1 Discrete Fourier Transform

For a set of complex numbers x_{n}, n = 0, …, N − 1, the DFT is defined as a set of complex numbers c_{k},

c_{k} = ∑_{n=0}^{N−1} x_{n}e^{−2πi \frac{nk}{N}} , k = 0, …, N − 1 .  (11.1)

The inverse DFT is given by

x_{n} = \frac{1}{N} ∑_{k=0}^{N−1} c_{k}e^{+2πi \frac{nk}{N}} .  (11.2)

These transformations can be viewed as expansion of the vector x_{n} in terms of the orthogonal basis of vectors e^{2πi \frac{kn}{N}},

∑_{n=0}^{N−1} (e^{2πi \frac{kn}{N}} ) (e^{−2πi \frac{k′n}{N}} ) = Nδ_{kk′} .  (11.3)

The DFT represent the amplitude and phase of the different sinusoidal components in the input data x_{n}.

The DFT is widely used in different fields, like spectral analysis, data compression, solution of partial differential equations and others.

11.1.1 Applications

Data compression

Several lossy (that is, with certain loss of data) image and sound compression methods employ DFT as an approximation for the Fourier series. The signal is discretized and transformed, and then the Fourier coefficients of high/low frequencies, which are assumed to be unnoticeable, are discarded. The decompressor computes the inverse transform based on this reduced number of Fourier coefficients.

Partial differential equations

Discrete Fourier transforms are often used to solve partial differential equations, where the DFT is used as an approximation for the Fourier series (which is recovered in the limit of infinite N). The advantage of this approach is that it expands the signal in complex exponentials e^{inx}, which are eigenfunctions of differentiation: \frac{d}{dx}e^{inx} = ine^{inx}. Thus, in the Fourier representation, differentiation is simply multiplication by in.

A linear differential equation with constant coefficients is transformed into an easily solvable algebraic equation. One then uses the inverse DFT to transform the result back into the ordinary spatial representation. Such an approach is called a spectral method.

Convolution and Deconvolution

FFT can be used to efficiently compute convolutions of two sequences. A convolution is the pairwise product of elements from two different sequences, such as in multiplying two polynomials or multiplying two long integers.

Another example comes from data acquisition processes where the detector introduces certain (typically Gaussian) blurring to the sampled signal. A reconstruction of the original signal can be obtained by deconvoluting the acquired signal with the detector’s blurring function.

11.2 Cooley-Tukey algorithm

In its simplest incarnation this algorithm re-expresses the DFT of size N = 2M in terms of two DFTs of size M,

c_{k} = ∑_{n=0}^{N−1} x_{n}e^{−2πi \frac{nk}{N}}

= ∑_{m=0}^{M−1} x_{2m}e^{−2πi \frac{mk}{M}} + e^{−2πi \frac{k}{N}} ∑_{m=0}^{M−1} x_{2m+1}e^{−2πi \frac{mk}{M}}

= \{ \begin{matrix}c_{k}^{(even)} + e^{−2πi \frac{k}{N}} c_{k}^{(odd)} & , k < M \\ c_{k−M}^{(even)} − e^{−2πi \frac{k−M}{N}} c_{k−M}^{(odd)} & , k ≥ M \end{matrix} ,  (11.4)

where c^{(even)} and c^{(odd)} are the DFTs of the even- and odd-numbered subsets of x.

This re-expression of a size-N DFT as two size-\frac{N}{2} DFTs is sometimes called the Danielson-Lanczos lemma. The exponents e^{−2πi \frac{k}{N}} are called twiddle factors.

The operation count by application of the lemma is reduced from the original N^{2} down to 2(N/2)^{2} + N/2 = N^{2}/2 + N/2 < N^{2}.

For N = 2^{p} Danielson-Lanczos lemma can be applied recursively until the data sets are reduced to one datum each. The number of operations is then reduced to O(N ln N) compared to the original O(N^{2}). The established library FFT routines, like FFTW and GSL, further reduce the operation count (by a constant factor) using advanced programming techniques like precomputing the twiddle factors, effective memory management and others.

11.3 Multidimensional DFT

For example, a two-dimensional set of data x_{n_{1}n_{2}}, n_{1} = 1 … N_{1}, n_{2} = 1 … N_{2} has the discrete Fourier transform

c_{k_{1}k_{2}} = ∑_{n_{1}=0}^{N_{1}−1} ∑_{n_{2}=0}^{N_{2}−1} x_{n_{1}n_{2}}e^{−2πi \frac{n_{1}k_{1}}{N_{1}}} e^{−2πi \frac{n_{2}k_{2}}{N_{2}}} .  (11.5)

11.4 C implementation

92 CHAPTER 11. FAST FOURIER TRANSFORM

Bibliography

[1] Hiroshi Akima. A new method of interpolation and smooth curve fitting based on local procedures. Journal of Assoc. for Comp. Mach., 17(4):589–602, 1970.

[2] Przemyslaw Bogacki and Lawrence F. Shampine. A 3(2) pair of Runge–Kutta formulas. Applied Mathematics Letters, 2(4):321–325, 1989.

[3] C.W. Clenshaw and A.R. Curtis. A method for numerical integration on an automatic computer. Numerische Mathematik, 2:197–205, 1960.

[4] M. Galassi et al. GNU Scientific Library Reference Manual. Network Theory Ltd, 3rd edition, 2009.

[5] Erwin Fehlberg. Low-order classical Runge-Kutta formulas with step size control and their application to some heat transfer problems. NASA Technical Report, 1969.

[6] Wallace Givens. Computation of plane unitary rotations transforming a general matrix to triangular form. J. SIAM, 6(1):26–50, 1958.

[7] Gene H. Golub. Some modified matrix eigenvalue problems. SIAM Rev., 15(2):318–334, 1973.

[8] Gene H. Golub and John H. Welsch. Calculation of Gauss quadrature rules. Mathematics of Computation, 23(106):221–230, 1969.

[9] Pedro Gonnet. Increasing the reliability of adaptive quadrature using explici interpolants. ACM Trans. Math. Soft., 37(3):26:2–26:32, 2010.

[10] A.S. Householder. Unitary triangularization of a nonsymmetric matrix. Journal of the ACM, 5(4):339–342, 1958.

[11] Aleksandr Semenovich Kronrod. Nodes and weights of quadrature formulas. Sixteen-place tables. Consultants Bureau, 1965.

[12] Masatake Mori. Quadrature formulas obtained by variable transformation and the DE–rule. Journal of Computational and Applied Mathematics, 12&13:119–130, 1985.

[13] T. N. L. Patterson. The optimum addition of points to quadrature formulae. Mathematics of Computation, 22(104):847–856, 1968.

[14] Robert Piessens, Elise de Doncker-Kapenga, Christoph W. Überhuber, and David Kahaner. QUADPACK: A subroutine package for automatic integration. Springer-Verlag, 1983.

[15] B. Riemann. Über die darstellbarkeit einer function durch eine trigonometrische reihe. Abhandlungen der Königlichen Gesellschaft der Wissenschaften zu Göttingen, 13:87–132, 1868.

Index

Akima sub-spline, 7 Arnoldi iteration, 86

back-substitution, 48 backtracking, 62 binary search, 3 Bogacki-Shampine method, 40 Butcher’s tableau, 39

characteristic equation, 81 classical quadrature, 16 cubature, 13 curve fitting, 1

Danielson-Lanczos lemma, 91 degree-n quadrature, 16

eigenvalue, 77 eigenvector, 77 embedded rule, 21 Euler’s method, 38

forward-substitution, 48

GMRES, 88

half-interval search, 3 Hessian matrix, 68 Householder reflection, 51 Householder transformation, 49, 51,

52

initial value problem, 35 interpolant, 1 inverse power iteration, 85

Jacobi matrix, 19 Jacobi operator, 19 Jacobi rotation, 78 Jacobian matrix, 62

Krylov matrix, 86 Krylov subspace, 86

Lagrange interpolating polynomial, 1 Lanczos iteration, 87 line search, 62 LU-decomposition, 54

matrix diagonalization, 77 mesh, 14

natural spline, 5 Newton’s step, 68 Newton-Cotes quadrature, 17 nodes, 13

orthogonal transformation, 77 overdetermined system, 57

partial pivoting, 55

QR-decomposition, 49 quadrature, 13

Rayleigh quotient, 86 rectangle rule, 14 Riemann integral, 13 Riemann sum, 13 root-finding, 61 rule, 13 Runge’s phenomenon, 1

secular equation, 81 similarity transformation, 77 spline, 2 96 INDEX trapezium rule, 14 triangular system, 48