The daily routine

I will present one simple algorithm per day. It will add up. In fact one can get by very nicely with about one dozen simple numerical recipes.

Day 1. Eulerian extrapolation

To solve a ﬁrst order differential equation

y ′(x) = f(x, y), y(0) = y0

we simply write it as

y(x--+--dx)----y(x)- = f(x, y(x)), y(x + dx) = y(x) + f (x,y(x)) ⋅ dx dx

which extrapolates the value of y(x + dx) from that of y(x). The increment dx should be small, and even then the algorithm is prone to accumulated errors.
To numerically solve

y ′′(x) = f (x, y(x), y′(x)), y(0) = y0, y′(0) = v0

we perform the Hamiltonian split, and get two ﬁrst order equations

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

and apply the ﬁrst method to both

y(x + dx) = y(x) + v(x) ⋅ dx, v(x + dx) = v(x) + f(x, y(x), v(x)) ⋅ dx

Day 2. Newton’s method

The technique is based on the equation

′ f (x) f (y) = f (x) + (y - x) f (x) + ⋅ ⋅⋅, y ≈ x - --′--- f (x)

The idea is that if x₀ is nearly a solution to f(x) = 0, then x₁ = x₀ - f(x0)-
f′(x0)

is even closer to solving the equation. This can be used to compute roots of functions, or to compute functions themselves to higher precision.
Example Let

X f (x) = ---- 1 x2

and the sequence

{x , x , x ,⋅⋅⋅}, x = x + f-(xn)-- 0 1 2 n+1 n f ′(xn)

converges on √ ---
X

.
The technique works for complex functions as well as real. Bisection may be superior for real-valued functions.
Warning Newton’s method can be unstable if x₀ is not chosen wisely.

Day 3. Numerical integration (naive)

A function f(x) can be numerically integrated over a domain a ≤ x ≤ b by coding the Riemann sum

∫ b (b - a) N∑ b - a f(x) dx = -------- ⋅ f (a + n ⋅ (------)) a N n=0 N

and by choosing N to be fairly large. This utterly fails for inﬁnite or semi-inﬁnite domains. Gaussian quadrature works in those cases, and is more efficient in all but the simplest cases.

Day 4. Naive sorting

Simple comparative sorting is slow, requiring n(n+1)
2 comparisons to sort n data, but it is sure-ﬁre. Heap sorting can be orders of magnitude faster, but it is more complex to encode. A well-coded heap sort program is somewhat inscrutible, but naive transposition sorting is straightforward.
To sort n values {x₀,x₁,x₂, ⋅⋅⋅ ,x_n-1} into ascending order of size, we run a pair of nested loops on the numbers stored in an array.

(for i=0;i<n-1;i++){
      for(j=i;j<n;j++){
           tmp=x[i];
           if(x[j]<x[i]){
                  x[i]=x[j];
                  x[j]=tmp;   /* we swap them */
           }
      }
}

Day 5. Numerical integration (stochastic)

Riemann sums are totally impractical for multi-dimensional integrals. For multiple integrals Monte Carlo or discrepancy sequencing may be the only practical algorithm for deﬁnite integration. Take the following steps to integrate a single-variable non-negative function f(x) over [a, b].
A. Pick a number Y _max > f(x) for all a ≤ x ≤ b.
B. Generate N pairs of random numbers (r_x,r_y) with a ≤ r_x ≤ b and 0 ≤ r_y ≤ Y _max.
C. for each pair, if r_y ≤ f(r_x), count the pair as a “hit”.
D. Let the number of “hits” be called h. Thinking of the rectangle of area Y _max ⋅ (b - a) as a target, the area ∫ _a^bf(x) dx as a “bull’s eye”, the ratio of the number of hits to the number of shots taken should be the area ratio, or

∫ b h f (x) dx = ---⋅ Ymax ⋅ (b - a) a N

Day 6. Special functions by recursion

Many special functions of mathematical physics possess three-term recursions, for example the Legendre polynomial P_n(x) for -1 ≤ x ≤ 1 obeys

(n + 1) Pn+1(x) = (2n + 1) Pn(x) - n Pn - 1(x), P0(x) = 1, P1(x) = x

Rather than writing a routine that calls itself recursively, which is very costly, a program to compute P_n(x) should return 1 if n = 0, x if n = 1, and if n > 1, the program should loop and build up successively higher-n values recursively.
A warning, many three-term recursions are numerically unstable to round-off errors performed in the machine. In this case one uses the Miller recursion to work down the chain to lower n values rather than higher.
Example Bessel functions obey

2n Jn+1(x) = ---Jn(x) - Jn - 1(x) x

which can be frightfully unstable. The Miller algorithm would work in the following way; to compute J₀(x), we begin high up, and pick

J6(x) = 0, J5(x) = 1

and go downwards;

J (x) = 10-J (x) - J (x) = 10- 4 x 5 6 x

J (x) = 8-J (x) - J (x) = 8-10-- 1 3 x 4 5 x x

6- 6- 8-10- 10- J2(x) = x J3(x) - J4(x) = x (x x - 1) - x

4 4 480 16 80 J1(x) = --J2(x) - J3(x) = --(------ ---) - ---+ 1 x x x3 x x2

and ﬁnally

2 2 1920 144 480 16 J0(x) = --J1(x) - J2(x) = --(------- -----+ 1) - -----+ --- x x x4 x2 x3 x

This result is not correct since in reality J₅(x)

1. Now we make use of the growth of the recursion; we can see that whatever x was, our value for J₀(x) is very large compared to the value for J₅(x) , and so it will clearly dominate in the well known sum

1 = J0(x) + 2J2(x) + 2J4(x) + 2J6(x) + ⋅⋅⋅

Each of our values for the various intervening Bessel functions is therefore off by a constant multiplicative factor

′′ ′′ 1 = J0(x) + 2J2(x) + 2J4(x) + 2J6(x) + ⋅⋅⋅

3840 768 18 480 16 10 = (------- -----+ ---) + 2( ------ ---) + 2( ---) + (0) x5 x3 x x3 x x

which we now divide by to restore our values to the correct values

3840-- 768 + 18 J0(x) ≈ ------------------x5------x3----x------------------ ( 38x450- 76x83 + 1x8) + 2( 48x03 - 1x6) + 2( 1x0) + (0)

For x = 1 this results in

3090-- J0(1) ≈ 4038 ≈ 0.76523

How close is this? The correct value is J₀(1) = 0.7651976866 ⋅⋅⋅

Day 7. Division-free matrix inversion

Newton’s method works for matrices as well as simple numerical types. Let

1 f (x) = (X - --) x

then

xn+1 = 2xn - xn ⋅ X ⋅ xn

and if x_n and X are matrices, the sequence {x₀, x₁, x₂, ⋅⋅⋅

} converges very rapidly on X^-1, and the entire computation is free of divisions. The matrix x₀ must be chosen very wisely. Deﬁne

∑ 2 ∣∣X ∣∣ = X i,j i,j

the sum of all matrix elements squared, and choose

x = --1---X T 0 ∣∣X ∣∣

and the sequence converges on X^-1 if X is invertible, or x₀ itself if X is non-invertible, faster than any other algorithm for medium to large matrices.

Day 8. Creating random doubles distributed with density f(x)

The accept-reject method of John Von Neumann is not the most efficient for every PDF, but it always works. To generate a random number n with PDF f(x);
A. Choose a domain, a ≤ x ≤ b,
B. for the PDF such f(x), and ﬁnd its largest value f_max on your domain.
C. Generate a uniformly distributed random number x on this domain. The code fragment below illustrates how to make this into a random ﬂoat between a and b.

#include <stdlib.h>
srand((unsigned long)getpid());
x=(b-a)*((double)rand()/(double)RAND_MAX)+a;

D. Pick a second uniformly distributed random number y on the interval 0 ≤ y ≤ f_max. If y ≤ f(x), then accept this number as being randomly distributed according to f(x). If y > f(x), reject it and pick a new pair (x,y). One look at the ﬁgure below tells you why this works; the probability of acceptance of x is the same as the probability that y < f(x), which is ℘(y) = ff(mxa)x .

Day 9. Least squares methods

Least squares methods are based on the use of minimax methods for expanding a sampled function f(x) into a sum basis functions {g_i(x)} that minimizes the standard deviation S;

M N ∑ ---1----∑ 2 f (x) = ai gi(x), S = N - 1 (yj - f (xj)) i=1 j=1

Example the nonlinear regression for the coefficients of the polynomial

y(x) = a + a x + a x2 0 1 2

that best ﬁts a set of data points

(x1, y1), (x2, y2),⋅⋅ ⋅,(xm, ym)

is gotten by minimizing the deviation

m 2 S = ---1----∑ (y - ∑ a xj )2 m - 1 k=1 k j=1 j k

with respect to the coefficients a₀,a₁ and a₂, resulting in the simultaneous linear equations

∑ ∑ ∑ 2 yk = na0 + a1 xk + a2 xk k k k

∑ ∑ ∑ 2 ∑ 3 ykxk = a0 xk + a1 x k + a2 xk k k k k

∑ ∑ ∑ ∑ ykx2k = a0 x2k + a1 x3k + a2 x4k k k k k

These can be solved for the coefficients a_i by writing these equations as a matrix relation and inverting the matrix.
Call S_x = ∑ _i=1^mx_i, S_x² = ∑ _i=1^mx_i², S_x³ = ∑ _i=1^mx_i³, S_y = ∑ _i=1^my_i, S_xy = ∑ _i=1^mx_i y_i and S_x²y = ∑ _i=1^mx_i² y_i and so on. We can write the equations above as

( ) ( ) ( ) | Sy | | n Sx Sx2 | | a0 | || Sxy || = || Sx Sx2 Sx3 || || a1 || ( ) ( ) ( ) Sx2y Sx2 Sx3 Sx4 a2

We can use our multiplicative inversion to solve numerically or

( ) ( ) - 1 ( ) | a0 | | n Sx Sx2 | | Sy | || a || = || S S 2 S 3 || || S || ( 1 ) ( x x x ) ( xy ) a2 Sx2 Sx3 Sx4 Sx2y

This method can be used to compute “Fourier transforms”, not true Fourier transforms, but numerically very good ﬁts to sampled functions. The technique is also the basis of the Gaussian quadrature method of integration.

Day 10. Runge-Kutta extrapolation

This is a technique for solving differential equations superior to Euler (or any other method for that matter).
Simple Runge-Kutta is a three-step extrapolation that works this way.
The equation y′′(x) = f(x,y)

1 yn+1 = yn + δx (y′ + --(k1 + 2k2)) n+1 6

y′ = y′ + 1(k + 4k + k ) n+1 n 6 1 2 3

in which

k1 = δx f (xn, yn)

δx- δx- ′ δx- k2 = δx f (xn + 2 ,yn + 2 yn + 8 k1)

k = δx f (x + δ , y + δx y′ + δx-k ) 3 n x n n 2 2

The equations y′(x) = f(x,y,z), z′(x) = g(x,y,z)

1- 1- yn+1 = yn + 2(k1 + k2), zn+1 = zn + 2( ℓ1 + ℓ2)

with

k1 = δx f(xn, yn, zn), ℓ1 = δx g(xn, yn,zn)

k2 = δx f (xn + δx, yn + k1, zn + ℓ1), ℓ2 = δx g(xn + δx, yn + k1,zn + ℓ1)

This is the scheme that is most useful in classical mechanics; we let x = t, y(x) = r(t) be the position versus time, and z(x) = p(t) be the momentum versus time.
This case is also quite easy to derive, and doing so illustrates the order of the error in the formula. We begin with the power series

2 ′ δ-- ′′ y(x + δ) = y(x) + δy (x) + 2 y (x) + ⋅⋅ ⋅

δ2 ∂ ′ ∂ ′ ∂ = y(x)+ δf (x, y,z)+ ---(---f (x,y, z) y (x)+ ---f (x,y, z) z (x)+ ---f (x, y,z))+ ⋅⋅⋅ 2 ∂y ∂z ∂x

and

2 z(x + δ) = z(x) + δz′(x) + δ--z′′(x) + ⋅⋅ ⋅ 2

δ2- ∂-- ′ ∂-- ′ = z(x) + δg(x, y, z) + ( g(x, y,z) y (x) + g(x, y,z) z (x) 2 ∂y ∂z

-∂-- + ∂x g(x, y,z)) + ⋅⋅ ⋅

We now just split off the second term and regroup

δ- δ- -∂- ′ ∂-- ′ y(x+ δ) = y(x)+ 2 f (x,y, z)+ 2[f (x, y,z)+ δ∂y f(x, y, z) y (x)+ δ ∂z f (x,y, z) z (x)

∂ + δ ---f (x, y,z) ] + ⋅⋅ ⋅ ∂x

Notice that the terms inside of the square brackets are in fact the ﬁrst few terms of the multi-variable Taylor series expansion for the function f(x,y,z)

∂-- ′ -∂- ′ -∂-- f (x, y,z) + δ f (x,y, z) y (x) + δ f (x, y,z) z (x) + δ f (x,y, z) + ⋅⋅⋅ ∂y ∂z ∂x

≈ f (x + δ,y + δy′,z + δz ′)

≈ f(x + δ, y + δf (x, y, z),z + δg(x, y, z))

and so we have

y(x + δ) = y(x) + δ(f (x, y,z) + f(x + δ,y + δf (x, y,z), z + δg(x, y, z))) + ⋅⋅ ⋅ 2

and similarly

δ- z(x + δ) = z(x) + (g(x, y, z) + g(x + δ,y + δf (x, y,z), z + δg(x, y, z))) + ⋅⋅ ⋅ 2

and so if we call

k1 = δf (x, y, z), k2 = δf (x + δ,y(x) + δf (x, y, z),z(x) + δg(x, y, z))

and

ℓ1 = δg(x, y, z), ℓ2 = δg(x + δ, y(x) + δf (x, y,z), z(x) + δg(x, y,z))

we have proven this Runge-Kutta formula. The others are done the same way, with more or less algebra.

Day 11. Eigenvalues of a matrix

Gershgorin’s theorem tells us that the magnitude of the eigenvalues of a real, symmetric matrix are less than the maximum value of the sum of absolute values of row elements. For example, given a matrix A, let

∑N r = ∣A ∣ i j=1 i,j

then if

ρ = M ax{r1, r2, ⋅⋅⋅, rN }

we have

∣λj ∣ ≤ ρ, ∀j

This number ρ is called the spectral radius of the matrix A.
Matrix eigenvalues can be extracted by the very crude but effective method of iteration using this fact. We construct an arbitrary vector v₀ and build up the sequence of vectors

vi = Avi - 1

which converges on the eigenvector of the largest eigenvalue, since if

Axi = λixi

is the equation satisﬁed by the i^th eigenvector, then we can expand

∑ v0 = aixi i

and we see that

∑ ∑ vk = ai(A)kxi = aiλkixi i i

and as k increases, this sum becomes dominated by the largest summand, corresponding to the largest magnitude eigenvalue. If for example the ℓ^th component of v_k is nonzero,

((A)k+1v0) ℓ λmax mag = lim -------k------ k→ ∞ ((A) v0)ℓ

A variation on this can be used to ﬁnd any eigenvalue, but it may be a great deal of work, and better methods do exist. Notice that if the eigenvalues of A are λ₁,λ₂, ⋅⋅⋅

,λ_N, then the eigenvalues of A - t ⋅ 1 are λ₁ - t,λ₂ - t, ⋅⋅⋅

,λ_N - t, and of (A - t ⋅ 1)^-1 are -1--
λ1-t

. If we choose t accordingly, the recursion

- 1 vi+1 = (A - t ⋅ 1) vi

will converge on the eigenvector of (A-t⋅ 1)^-1 of largest magnitude, giving the eigenvalue of A closest to the number t. If A has close eigenvalues, this method fails.

The advantages of this algorithm is that it is easy to understand and to encode. We need only a matrix multipliction routine, and one to rescale the matrix. Other more powerful methods are bigger coding jobs and involve more complex linear algebra principles.