Solutions of linear equations

Conjugate gradient method

The conjugate gradient (CG) method solves $Ax = b$ for a symmetric, positive definite matrix $A$. It is equivalent to minimising the quadratic energy

$$\Pi(x) = \tfrac{1}{2}x^T A x - b^T x$$

whose minimum is exactly the solution of $Ax=b$. The minimum is searched along directions $d_k$ that are A-conjugate, i.e. $d_i^T A d_j = 0$ for $i \ne j$. Each step minimises $\Pi$ along one direction without spoiling the earlier ones, so in exact arithmetic the solution is reached in at most $n$ steps for an $n \times n$ system. The residual $r = b - Ax$ is the negative gradient of $\Pi$.

Compared with direct elimination, CG needs only matrix-vector products, so it stores only the non-zero terms of $A$. It is good for large sparse systems, and rounding errors do not accumulate in a long elimination.

Algorithm

1. Choose a starting vector $x_0$ (usually $0$). Compute the residual $r_0 = b - Ax_0$ and set the first search direction $d_0 = r_0$.
2. For $k = 0, 1, 2, \dots$ repeat:
   • Step length: $\alpha_k = \dfrac{r_k^T r_k}{d_k^T A d_k}$
   • Update the solution: $x_{k+1} = x_k + \alpha_k d_k$
   • Update the residual: $r_{k+1} = r_k - \alpha_k A d_k$
   • Check convergence: stop if $\|r_{k+1}\| <$ tolerance.
   • Direction factor: $\beta_k = \dfrac{r_{k+1}^T r_{k+1}}{r_k^T r_k}$
   • New direction: $d_{k+1} = r_{k+1} + \beta_k d_k$
3. The last $x_{k+1}$ is the solution.

Convergence is faster if $A$ has clustered eigenvalues, and preconditioning is used to improve this.

Matrix form: $A=\begin{bmatrix}2&-1&0\\-1&2&-1\\0&-1&2\end{bmatrix}$, $b=\{4,0,0\}^T$. $A$ is symmetric and positive definite. Start with $x_0=\{0,0,0\}^T$.

$r_0=b-Ax_0=\{4,0,0\}^T$, $d_0=r_0$.

Iteration 1

$$\begin{aligned} Ad_0&=\{8,-4,0\}^T\\ \alpha_0&=\frac{r_0^Tr_0}{d_0^TAd_0}=\frac{16}{32}=0.5\\ x_1&=\{2,0,0\}^T\\ r_1&=r_0-\alpha_0Ad_0=\{0,2,0\}^T\\ \beta_0&=\frac{4}{16}=0.25,\quad d_1=r_1+\beta_0d_0=\{1,2,0\}^T \end{aligned}$$

Iteration 2

$$\begin{aligned} Ad_1&=\{0,3,-2\}^T\\ \alpha_1&=\frac{4}{6}=0.6667\\ x_2&=\{2.6667,\,1.3333,\,0\}^T\\ r_2&=\{0,0,1.3333\}^T\\ \beta_1&=\frac{1.7778}{4}=0.4444,\quad d_2=\{0.4444,\,0.8889,\,1.3333\}^T \end{aligned}$$

Iteration 3

$$\begin{aligned} Ad_2&=\{0,0,1.7778\}^T\\ \alpha_2&=\frac{1.7778}{2.3704}=0.75\\ x_3&=\{3,\,2,\,1\}^T\\ r_3&=\{0,0,0\}^T \end{aligned}$$

The residual is zero after $n=3$ iterations, as expected for CG.

Answer: $x_1=3,\ x_2=2,\ x_3=1$. (Check: $2(3)-2=4$, $-3+4-1=0$, $-2+2=0$.)

$A=\begin{bmatrix}3&0&2\\0&1&1\\2&1&3\end{bmatrix}$ is symmetric and positive definite (leading minors $3,\ 3,\ 4>0$). $b=\{1,0,-1\}^T$. Take $x_0=0$: $r_0=b$, $d_0=r_0=\{1,0,-1\}^T$.

Iteration 1

$$\begin{aligned} Ad_0&=\{1,\,-1,\,-1\}^T\\ \alpha_0&=\frac{r_0^Tr_0}{d_0^TAd_0}=\frac{2}{2}=1\\ x_1&=\{1,\,0,\,-1\}^T\\ r_1&=r_0-\alpha_0Ad_0=\{0,\,1,\,0\}^T\\ \beta_0&=\frac{1}{2}=0.5,\quad d_1=r_1+\beta_0d_0=\{0.5,\,1,\,-0.5\}^T \end{aligned}$$

Iteration 2

$$\begin{aligned} Ad_1&=\{0.5,\,0.5,\,0.5\}^T\\ \alpha_1&=\frac{1}{0.5}=2\\ x_2&=\{2,\,2,\,-2\}^T\\ r_2&=\{-1,\,0,\,-1\}^T\\ \beta_1&=\frac{2}{1}=2,\quad d_2=\{0,\,2,\,-2\}^T \end{aligned}$$

Iteration 3

$$\begin{aligned} Ad_2&=\{-4,\,0,\,-4\}^T\\ \alpha_2&=\frac{2}{8}=0.25\\ x_3&=\{2,\,2.5,\,-2.5\}^T\\ r_3&=\{0,\,0,\,0\}^T \end{aligned}$$

Answer: $x_1=2,\ x_2=2.5,\ x_3=-2.5$. Check: $3(2)+2(-2.5)=1$, $2.5-2.5=0$, $2(2)+2.5-7.5=-1$.

A system $[A]\{x\}=\{b\}$ is solved by direct or iterative methods.

Direct methods

They give the exact answer (apart from round-off) after a fixed number of operations.

• Gauss elimination: reduce $A$ to upper triangular form, then back-substitute.
• LU decomposition / Cholesky ($A=LL^T$ for symmetric positive definite matrices): factor once, then solve for several right-hand sides.
• Gauss-Jordan and matrix inverse.

Example: $2x+y=5,\ x+3y=10$. Eliminating, $R_2\to R_2-\tfrac12R_1$ gives $2.5y=7.5$, so $y=3$, $x=1$.

Iterative methods

Start from a guess and improve it until the change is small. They suit large sparse systems.

• Jacobi: $x_i^{(k+1)}=\dfrac{1}{a_{ii}}\Big(b_i-\sum_{j\ne i}a_{ij}x_j^{(k)}\Big)$
• Gauss-Seidel: same but uses new values immediately.
• SOR: Gauss-Seidel with relaxation factor $\omega$.
• Conjugate gradient: minimises $\tfrac12x^TAx-b^Tx$ for symmetric positive definite $A$.

Example (Gauss-Seidel): $4x+y=9,\ x+3y=7$, start $(0,0)$.

• Iteration 1: $x=9/4=2.25$, $y=(7-2.25)/3=1.583$.
• Iteration 2: $x=(9-1.583)/4=1.854$, $y=(7-1.854)/3=1.715$.

The values move to the exact solution $x=20/11=1.818$, $y=19/11=1.727$.

Comparison

Sparse matrix: a matrix in which most entries are zero. Stiffness matrices of FEM are sparse because a node is connected only to the nodes of neighbouring elements.

Banded matrix: a matrix whose non-zero entries lie in a narrow band about the main diagonal, i.e. $a_{ij}=0$ for $|i-j|>b$. The half-bandwidth is

$$b=(\text{max difference of node numbers in an element}+1)\times \text{dof per node}$$

A symmetric matrix of order $n$ and half-bandwidth $b$ is stored with only $n\times b$ numbers instead of $n^2$.

 x x . . .
 x x x . .      band along the
 . x x x .      diagonal, zeros
 . . x x x      outside
 . . . x x

Memory optimisation: reducing the memory needed to store and solve $[K]\{U\}=\{F\}$. Techniques are:

• store only the upper triangle (symmetry);
• store only the band (banded storage) or the skyline profile;
• store only non-zero terms with index arrays (sparse storage);
• renumber nodes to minimise bandwidth. Less memory also gives fewer operations: the cost of banded elimination is about $n b^2$ instead of $n^3$.

Iterative methods start from an assumed solution and improve it repeatedly until the change is within tolerance. They are used for large sparse systems $[A]\{x\}=\{b\}$ and need diagonal dominance (or symmetry and positive definiteness) for convergence.

1. Jacobi method: every unknown is computed using values of the previous iteration only.

$$x_i^{(k+1)}=\frac{1}{a_{ii}}\Big(b_i-\sum_{j\ne i}a_{ij}x_j^{(k)}\Big)$$

2. Gauss-Seidel method: newly computed values are used immediately in the same iteration:

$$x_i^{(k+1)}=\frac{1}{a_{ii}}\Big(b_i-\sum_{j<i}a_{ij}x_j^{(k+1)}-\sum_{j>i}a_{ij}x_j^{(k)}\Big)$$

It converges about twice as fast as Jacobi. 3. Successive over-relaxation (SOR): $x_i^{new}=(1-\omega)x_i^{old}+\omega\,x_i^{GS}$ with $1<\omega<2$ to speed up convergence. 4. Conjugate gradient method: for symmetric positive definite matrices; minimises the quadratic energy along conjugate directions and converges in at most $n$ steps.

Direct and iterative methods

• Direct methods give the solution in a finite, known number of operations (exact apart from round-off). Examples: Gauss elimination, Gauss-Jordan, LU and Cholesky decomposition, matrix inverse (Cramer's rule). They are best for small or dense systems.
• Iterative methods start from a trial vector and improve it repeatedly until the change is below a tolerance. Examples: Jacobi, Gauss-Seidel, SOR, conjugate gradient. They are best for large sparse systems and need convergence conditions such as diagonal dominance.

Gauss-Seidel solution

The system is diagonally dominant ($12>8$, $5>4$, $13>10$), so it converges. Rearranging:

$$x=\frac{1-3y+5z}{12},\quad y=\frac{28-x-3z}{5},\quad z=\frac{76-3x-7y}{13}$$

Start $(x,y,z)=(1,2,3)$.

Example of iteration 1: $x=(1-3(2)+5(3))/12=0.8333$; $y=(28-0.8333-9)/5=3.6333$; $z=(76-2.5-25.4333)/13=3.6974$.

Answer: $x\approx1,\ y\approx3,\ z\approx4$ (exact solution $x=1,y=3,z=4$).

In FEM the stiffness matrix $[K]$ is symmetric, sparse and has non-zero terms only near the diagonal. Storing the full $n\times n$ matrix wastes memory.

Banded storage

Only the band is stored: the diagonal and the $b-1$ upper terms of each row (using symmetry). Half-bandwidth $b=(\text{max node difference in an element}+1)\times$ dof per node.

Example: a bar of 100 nodes (1 dof), elements connect nodes $i$ and $i+1$, so $b=2$. Full storage $=100^2=10000$ numbers; banded storage $=100\times2=200$ numbers.

 Full 5x5 (symmetric)        Banded (b=2)
 a11 a12  0   0   0          a11 a12
 a12 a22 a23  0   0          a22 a23
  0  a23 a33 a34  0   ->     a33 a34
  0   0  a34 a44 a45         a44 a45
  0   0   0  a45 a55         a55  *

Skyline (variable band / profile) storage

Each column is stored from the first non-zero entry down to the diagonal, so the band height varies with column. A pointer array stores the diagonal positions.

Example: column heights 1, 2, 2, 3, 1 give 9 stored terms against 15 for the full upper triangle of a $5\times5$ matrix.

Skyline storage stores fewer terms than a constant band when a few nodes have a large connectivity. Fill-in during elimination remains inside the skyline, so the factorisation also needs no extra storage. Node renumbering (e.g. reverse Cuthill-McKee) reduces the bandwidth and so the memory and operation count.

Algorithm

1. Choose a starting vector $x_0$ (usually $0$). Compute the residual $r_0 = b - Ax_0$ and set the first search direction $d_0 = r_0$.
2. For $k = 0, 1, 2, \dots$ repeat:
   • Step length: $\alpha_k = \dfrac{r_k^T r_k}{d_k^T A d_k}$
   • Update the solution: $x_{k+1} = x_k + \alpha_k d_k$
   • Update the residual: $r_{k+1} = r_k - \alpha_k A d_k$
   • Check convergence: stop if $\|r_{k+1}\| <$ tolerance.
   • Direction factor: $\beta_k = \dfrac{r_{k+1}^T r_{k+1}}{r_k^T r_k}$
   • New direction: $d_{k+1} = r_{k+1} + \beta_k d_k$
3. The last $x_{k+1}$ is the solution.

Limitations

• Applies only to symmetric, positive definite matrices (not to general or indefinite systems without modification).
• Convergence slows for ill-conditioned matrices (large ratio of largest to smallest eigenvalue), so preconditioning is often needed.
• Round-off errors destroy the conjugacy of the directions, so more than $n$ iterations may be required in practice.
• Gives one right-hand side at a time; a direct (factored) solver is more efficient for many load cases.
• Needs a good stopping criterion; residual may fall while the error is still large.

Linear system $[A]\{x\}=\{b\}$ is solved by:

1. Direct methods (finite number of steps, exact apart from round-off): Gauss elimination (reduce to upper triangular form, back-substitute), Gauss-Jordan, LU decomposition, Cholesky decomposition $A=LL^T$ for symmetric positive definite $A$, and the matrix inverse. In FEM they use banded or skyline storage.
2. Iterative methods (start from a guess and improve): Jacobi, Gauss-Seidel, successive over-relaxation, and the conjugate gradient method. They suit very large sparse systems since the matrix is never changed and no fill-in occurs. Jacobi and Gauss-Seidel need diagonal dominance for convergence; CG needs a symmetric positive definite matrix.

Conjugate gradient, two iterations

$A=\begin{bmatrix}12&-6&0\\-6&12&-6\\0&-6&6\end{bmatrix}$, $b=\{24,24,0\}^T$, $x_0=\{4,4,0\}^T$.

$Ax_0=\{12(4)-6(4),\,-6(4)+12(4),\,-6(4)\}=\{24,24,-24\}^T$, so

$$r_0=b-Ax_0=\{0,0,24\}^T,\quad d_0=r_0$$

Iteration 1

$$\begin{aligned} Ad_0&=\{0,-144,144\}^T\\ \alpha_0&=\frac{r_0^Tr_0}{d_0^TAd_0}=\frac{576}{3456}=0.1667\\ x_1&=x_0+\alpha_0d_0=\{4,\,4,\,4\}^T\\ r_1&=r_0-\alpha_0Ad_0=\{0,\,24,\,0\}^T\\ \beta_0&=\frac{576}{576}=1,\quad d_1=r_1+\beta_0d_0=\{0,\,24,\,24\}^T \end{aligned}$$

Iteration 2

$$\begin{aligned} Ad_1&=\{-144,\,144,\,0\}^T\\ \alpha_1&=\frac{576}{3456}=0.1667\\ x_2&=x_1+\alpha_1d_1=\{4,\,8,\,8\}^T\\ r_2&=r_1-\alpha_1Ad_1=\{24,\,0,\,0\}^T \end{aligned}$$

Answer: after two iterations $x^{(2)}=\{4,\,8,\,8\}^T$ with residual $\{24,0,0\}^T$. A third iteration gives the exact solution $x=\{8,12,12\}^T$ (the system has $n=3$ unknowns).

Discrete Fourier Transform (DFT)

The DFT converts $N$ equally spaced samples $x_n$ of a signal in the time domain into $N$ complex frequency components $X_k$:

$$X_k=\sum_{n=0}^{N-1}x_n\,e^{-i2\pi kn/N},\quad k=0,\dots,N-1$$

and the inverse is

$$x_n=\frac{1}{N}\sum_{k=0}^{N-1}X_k\,e^{i2\pi kn/N}$$

$|X_k|$ gives the amplitude and $\arg X_k$ the phase of the component at frequency $f_k=k/(N\Delta t)$. Direct evaluation needs $N^2$ complex multiplications.

Fast Fourier Transform (FFT)

The FFT (Cooley and Tukey, 1965) is an efficient algorithm for computing the DFT. It splits the sum into even- and odd-numbered samples repeatedly (for $N=2^m$):

$$X_k=E_k+W_N^{k}\,O_k,\qquad X_{k+N/2}=E_k-W_N^{k}\,O_k,\qquad W_N=e^{-i2\pi/N}$$

This "butterfly" operation reduces the work to $\tfrac{N}{2}\log_2N$ multiplications, so $N=1024$ needs about 5,000 operations instead of 10^6.

Uses in civil engineering: frequency content of earthquake records, vibration of bridges and buildings, spectral analysis of waves, signal filtering.

$A=\begin{bmatrix}10&-6&0\\-6&8&-2\\0&-2&5\end{bmatrix}$ (symmetric, positive definite), $b=\{12,0,0\}^T$. Take $x_0=\{0,0,0\}^T$, so $r_0=b=\{12,0,0\}^T$, $d_0=r_0$.

Iteration 1

$$\begin{aligned} Ad_0&=\{120,-72,0\}^T\\ \alpha_0&=\frac{144}{1440}=0.1\\ x_1&=\{1.2,\,0,\,0\}^T\\ r_1&=\{0,\,7.2,\,0\}^T\\ \beta_0&=\frac{51.84}{144}=0.36,\quad d_1=\{4.32,\,7.2,\,0\}^T \end{aligned}$$

Iteration 2

$$\begin{aligned} Ad_1&=\{0,\,31.68,\,-14.4\}^T\\ \alpha_1&=\frac{51.84}{228.10}=0.22727\\ x_2&=\{2.1818,\,1.6364,\,0\}^T\\ r_2&=\{0,\,0,\,3.2727\}^T\\ \beta_1&=\frac{10.7107}{51.84}=0.20661,\quad d_2=\{0.8926,\,1.4876,\,3.2727\}^T \end{aligned}$$

Iteration 3

$$\begin{aligned} Ad_2&=\{0,\,0,\,13.388\}^T\\ \alpha_2&=\frac{10.7107}{43.815}=0.24444\\ x_3&=\{2.4,\,2.0,\,0.8\}^T\\ r_3&=\{0,\,0,\,0\}^T \end{aligned}$$

Answer: $x_1=2.4,\ x_2=2.0,\ x_3=0.8$. Check: $10(2.4)-6(2)=12$, $-6(2.4)+8(2)-2(0.8)=0$, $-2(2)+5(0.8)=0$.

Algorithm for the conjugate gradient method

1. Choose a starting vector $x_0$ (usually $0$). Compute the residual $r_0 = b - Ax_0$ and set the first search direction $d_0 = r_0$.
2. For $k = 0, 1, 2, \dots$ repeat:
   • Step length: $\alpha_k = \dfrac{r_k^T r_k}{d_k^T A d_k}$
   • Update the solution: $x_{k+1} = x_k + \alpha_k d_k$
   • Update the residual: $r_{k+1} = r_k - \alpha_k A d_k$
   • Check convergence: stop if $\|r_{k+1}\| <$ tolerance.
   • Direction factor: $\beta_k = \dfrac{r_{k+1}^T r_{k+1}}{r_k^T r_k}$
   • New direction: $d_{k+1} = r_{k+1} + \beta_k d_k$
3. The last $x_{k+1}$ is the solution.

The method needs a symmetric positive definite matrix and gives the exact solution in at most $n$ steps.

Gauss-Seidel solution

The matrix is diagonally dominant ($2>1$, $6>3$, $8>7$), so Gauss-Seidel converges. Rewrite the equations:

$$x_1=\frac{2+x_2}{2},\qquad x_2=\frac{-4-x_1+2x_3}{6},\qquad x_3=\frac{5-4x_1+3x_2}{8}$$

Start with $x^{(0)}=(0,0,0)$ and use the newest values at once.

Example of iteration 1: $x_1=(2+0)/2=1$; $x_2=(-4-1+0)/6=-0.8333$; $x_3=(5-4+3(-0.8333))/8=-0.1875$.

Answer: $x_1\approx0.62,\ x_2\approx-0.76,\ x_3\approx0.03$.

Linear system $[A]\{x\}=\{b\}$ is solved by:

1. Direct methods (finite number of steps, exact apart from round-off): Gauss elimination (reduce to upper triangular form, back-substitute), Gauss-Jordan, LU decomposition, Cholesky decomposition $A=LL^T$ for symmetric positive definite $A$, and the matrix inverse. In FEM they use banded or skyline storage.
2. Iterative methods (start from a guess and improve): Jacobi, Gauss-Seidel, successive over-relaxation, and the conjugate gradient method. They suit very large sparse systems since the matrix is never changed and no fill-in occurs. Jacobi and Gauss-Seidel need diagonal dominance for convergence; CG needs a symmetric positive definite matrix.

Algorithm of the conjugate gradient method

The conjugate gradient (CG) method solves $Ax = b$ for a symmetric, positive definite matrix $A$. It is equivalent to minimising the quadratic energy

$$\Pi(x) = \tfrac{1}{2}x^T A x - b^T x$$

whose minimum is exactly the solution of $Ax=b$. The minimum is searched along directions $d_k$ that are A-conjugate, i.e. $d_i^T A d_j = 0$ for $i \ne j$. Each step minimises $\Pi$ along one direction without spoiling the earlier ones, so in exact arithmetic the solution is reached in at most $n$ steps for an $n \times n$ system. The residual $r = b - Ax$ is the negative gradient of $\Pi$.

1. Choose a starting vector $x_0$ (usually $0$). Compute the residual $r_0 = b - Ax_0$ and set the first search direction $d_0 = r_0$.
2. For $k = 0, 1, 2, \dots$ repeat:
   • Step length: $\alpha_k = \dfrac{r_k^T r_k}{d_k^T A d_k}$
   • Update the solution: $x_{k+1} = x_k + \alpha_k d_k$
   • Update the residual: $r_{k+1} = r_k - \alpha_k A d_k$
   • Check convergence: stop if $\|r_{k+1}\| <$ tolerance.
   • Direction factor: $\beta_k = \dfrac{r_{k+1}^T r_{k+1}}{r_k^T r_k}$
   • New direction: $d_{k+1} = r_{k+1} + \beta_k d_k$
3. The last $x_{k+1}$ is the solution.

Gaussian elimination is a direct method that reduces $[A]$ to triangular form, while the conjugate gradient (CG) method is an iterative method that minimises $\tfrac12x^TAx-b^Tx$ along A-conjugate directions. CG is preferred for large problems for these reasons:

1. Memory: CG needs only matrix-vector products $Ad$, so only the non-zero terms of the sparse matrix are stored. Gaussian elimination produces fill-in, which destroys sparsity inside the band.
2. Operation count: each CG iteration costs about $n\times(\text{non-zeros per row})$ operations, and a good solution often comes in far fewer than $n$ iterations, whereas elimination costs about $n^3/3$ (or $nb^2$ for a banded matrix).
3. Accuracy: round-off errors in elimination accumulate over $n$ steps; in CG each iteration uses the original matrix and corrects the error.
4. Flexibility: the solution can be stopped when the residual is small enough, and a good starting guess can be used; elimination cannot be stopped early.
5. Parallelism: matrix-vector products are easy to parallelise, and the matrix need not even be assembled (element-by-element products).
6. In exact arithmetic CG finishes in at most $n$ steps and with preconditioning converges faster.

Limitation: CG is applicable only to symmetric positive definite matrices, and for several load vectors a factorised direct solver may be cheaper.

$A=\begin{bmatrix}3&0&1\\0&-1&3\\1&3&0\end{bmatrix}$, $b=\{1,-12,2\}^T$. $A$ is symmetric (but not positive definite, as its eigenvalues are $3.61,\,2,\,-3.61$). The CG recurrences can still be applied as long as no denominator $d^TAd$ is zero; $\alpha$ may then become negative. Take $x_0=0$, so $r_0=b$ and $d_0=r_0=\{1,-12,2\}^T$.

Iteration 1

$$\begin{aligned} Ad_0&=\{5,\,18,\,-35\}^T\\ \alpha_0&=\frac{149}{-281}=-0.53025\\ x_1&=\{-0.5302,\,6.3630,\,-1.0605\}^T\\ r_1&=\{3.6512,\,-2.4555,\,-16.5587\}^T\\ \beta_0&=\frac{293.52}{149}=1.97015,\quad d_1=\{5.6214,\,-26.0973,\,-12.6184\}^T \end{aligned}$$

Iteration 2

$$\begin{aligned} Ad_1&=\{4.2458,\,-11.7579,\,-72.6706\}^T\\ \alpha_1&=\frac{293.52}{1247.5}=0.23527\\ x_2&=\{0.7923,\,0.2230,\,-4.0293\}^T\\ r_2&=\{2.6523,\,0.3108,\,0.5388\}^T\\ \beta_1&=0.025282,\quad d_2=\{2.7944,\,-0.3490,\,0.2197\}^T \end{aligned}$$

Iteration 3

$$\begin{aligned} Ad_2&=\{8.6031,\,1.0082,\,1.7475\}^T\\ \alpha_2&=0.30830\\ x_3&=\{1.6538,\,0.1154,\,-3.9615\}^T\\ r_3&\approx\{0,\,0,\,0\} \end{aligned}$$

The residual is zero after three iterations, so iterations 4 and 5 are not needed.

Answer: $x_1=1.6538\ (43/26),\ x_2=0.1154\ (3/26),\ x_3=-3.9615\ (-103/26)$. Check: $3(1.6538)-3.9615=1.0$, $-0.1154-11.8846=-12$, $1.6538+0.3462=2$.

Fast Fourier transform: idea

The DFT of $N$ samples is $X_k=\sum_{n=0}^{N-1}x_n W_N^{kn}$ with $W_N=e^{-i2\pi/N}$, needing $N^2$ operations. The Cooley-Tukey radix-2 FFT (for $N=2^m$) splits the series into even and odd samples:

$$X_k=E_k+W_N^kO_k,\qquad X_{k+N/2}=E_k-W_N^kO_k,\quad k=0,\dots,\tfrac N2-1$$

and repeats the split $\log_2N$ times. The cost is $\tfrac N2\log_2N$ complex multiplications.

Algorithm

1. Read $N$ (a power of 2) and the samples $x_n$ at interval $\Delta t$.
2. Reorder the samples in bit-reversed order.
3. For stage length $len=2,4,8,\dots,N$: compute $W=e^{-i2\pi/len}$; for each block, for each $k=0\dots len/2-1$ form $u=x_{i+k}$, $v=x_{i+k+len/2}\,W^k$, then set $x_{i+k}=u+v$, $x_{i+k+len/2}=u-v$.
4. After the last stage $x$ holds $X_k$. Amplitude $=2|X_k|/N$, frequency $f_k=k/(N\Delta t)$.

 x0 --o----o----o-- X0
 x4 --o----|    |
 x2 --o----o    |   (butterflies:
 x6 --o----|    |    a+Wb , a-Wb)
 ...

Program (C)

#include <stdio.h>
#include <math.h>
#include <complex.h>
#define PI 3.14159265358979323846

/* in-place radix-2 FFT, n must be a power of 2 */
void fft(double complex x[], int n)
{
    for (int i = 1, j = 0; i < n; i++) {      /* bit reversal */
        int bit = n >> 1;
        for (; j & bit; bit >>= 1) j ^= bit;
        j ^= bit;
        if (i < j) { double complex t = x[i]; x[i] = x[j]; x[j] = t; }
    }
    for (int len = 2; len <= n; len <<= 1) {  /* butterfly stages */
        double complex wlen = cexp(-2.0 * I * PI / len);
        for (int i = 0; i < n; i += len) {
            double complex w = 1.0;
            for (int k = 0; k < len / 2; k++) {
                double complex u = x[i + k];
                double complex v = x[i + k + len / 2] * w;
                x[i + k]           = u + v;
                x[i + k + len / 2] = u - v;
                w *= wlen;
            }
        }
    }
}

int main(void)
{
    int n = 8; double fs = 8.0; double complex x[8];
    for (int i = 0; i < n; i++) {             /* 1 Hz (amp 2) + 3 Hz (amp 1) */
        double t = i / fs;
        x[i] = 2*sin(2*PI*1*t) + 1*sin(2*PI*3*t);
    }
    fft(x, n);
    for (int k = 0; k < n; k++)
        printf("f=%.1f Hz  amp=%.3f\n", k*fs/n, 2*cabs(x[k])/n);
    return 0;
}

Example: time domain and frequency domain

Signal $x(t)=2\sin(2\pi t)+\sin(6\pi t)$ sampled at 8 Hz for 1 s ($N=8$). The program output is:

(components at $k=5,6,7$ are the mirror image of $k=3,2,1$.)

• Time domain shows the amplitude of vibration against time, peak value, duration, arrival time and decay of the signal, but the individual components are mixed together.
• Frequency domain shows the dominant frequencies (1 Hz and 3 Hz), their amplitudes (2 and 1) and phase. For a structure, the peak frequency is the natural frequency, so resonance with earthquake or machine frequency can be checked; noise is separated by filtering.
• Frequency resolution is $\Delta f=1/(N\Delta t)=1$ Hz; the highest frequency identified is the Nyquist frequency $f_s/2=4$ Hz.