CHAPTER 2: SOLUTION OF LINEAR EQUATIONS

2.2 SOLUTION OF SYSTEM OF LINEAR EQUATIONS

In computational engineering (such as Finite Element Analysis), the physical problem ultimately reduces to solving a large system of simultaneous algebraic linear equations of the form:

Where \([A]\) is the known square coefficient matrix (e.g., stiffness matrix), \(\{b\}\) is the known vector of constants (e.g., applied forces), and \(\{x\}\) is the unknown vector to be solved (e.g., nodal displacements). The solution techniques for these equations are broadly classified into two categories: Direct Methods and Iterative Methods.

1. Direct Methods

Direct methods attempt to solve the system of equations by a finite sequence of standard algebraic operations. In the absence of mathematical round-off errors, a direct method will yield the exact analytical solution.

• Gauss Elimination: Reduces the matrix \([A]\) to an upper triangular matrix \([U]\) using forward elimination, followed by backward substitution to find \(\{x\}\).
• LU Decomposition: Factors the matrix \([A]\) into the product of a lower triangular matrix \([L]\) and an upper triangular matrix \([U]\) such that \([A] = [L][U]\). It is highly efficient when solving for multiple different load vectors \(\{b\}\).
• Cholesky Factorization: A specialized decomposition \([A] = [L][L]^T\) used exclusively when \([A]\) is symmetric and positive-definite (common in structural stiffness matrices).

2. Iterative Methods

Iterative methods start with an initial "guess" for the solution vector \(\{x\}^{(0)}\) and successively refine this guess through a looped algorithm until the solution converges within an acceptable tolerance margin (\( \epsilon \)).

• Jacobi Method: Solves each equation for one specific variable assuming all other variables take their values from the previous iteration step.
• Gauss-Seidel Method: Similar to Jacobi, but uses the most recently updated values of variables immediately within the current iteration step, typically leading to faster convergence.
• Relaxation Methods: Introduces an extrapolation factor (like Successive Over-Relaxation, SOR) to artificially accelerate the convergence rate of the Gauss-Seidel method.

2.3 DATA STORAGE AND MEMORY OPTIMIZATION

Sparse and Banded Matrices

In structural engineering, the resulting global stiffness matrix \([K]\) for a large structure typically contains thousands of rows and columns.

• Sparse Matrix: A matrix populated primarily with zeros. Since each node in a structure is only connected to a few adjacent nodes, the vast majority of entries in the stiffness matrix are zero.
• Banded Matrix: A sparse matrix where all non-zero elements are concentrated within a narrow diagonal band around the main diagonal. The width of this non-zero band is dictated by the node numbering sequence of the mesh.

Memory Optimization Techniques

Standard computer storage of an \(n \times n\) matrix requires allocating memory for \(n^2\) elements. For a 10,000-degree-of-freedom system, this requires storing 100,000,000 numbers, which overwhelms RAM, even though 95% of those numbers are just zero. Memory optimization strategies aim to store only the non-zero elements.

1. Banded Matrix Storage

Instead of storing the full \(n \times n\) square, we only store a rectangular array of dimension \(n \times m_B\), where \(m_B\) is the half-bandwidth. All zeroes outside this diagonal band are completely ignored by the computer's memory.

2. Skyline Storage Scheme (Variable Bandwidth)

Even within the diagonal band, there may be many zeros. The skyline storage scheme improves upon banded storage. If there are zeroes at the top of a column, operations and storage are only required for the elements starting from the first non-zero value down to the main diagonal.

The jagged boundary separating the top zeroes from the first non-zero element in each column is visually termed the skyline. Only the numerical values "below the skyline" are stored in a continuous 1D vector array, dramatically minimizing memory consumption.

2.4 CONJUGATE GRADIENT METHOD

The Conjugate Gradient Method (CGM) is an advanced mathematical algorithm for the numerical solution of particular systems of linear equations, specifically those whose matrix is symmetric and positive-definite. It is an iterative method that finds the exact solution of an \(n \times n\) matrix in at most \(n\) steps, but is typically stopped much earlier when a desired tolerance is reached.

Algorithm for Conjugate Gradient Method

To solve the system \( [A]\{x\} = \{b\} \):

1. Initialization:
   • Choose an initial guess vector \( \{x\}_{0} \) (often just a vector of zeros).
   • Compute the initial residual vector: \( \{r\}_{0} = \{b\} - [A]\{x\}_{0} \).
   • Set the initial search direction vector equal to the residual: \( \{d\}_{0} = \{r\}_{0} \).
2. Iterative Loop: For \( k = 0, 1, 2, \dots \) until convergence:
   1. Compute the optimal step size \( \alpha_k \) along the search direction:
      $$ \alpha_k = \frac{\{r\}_k^T \{r\}_k}{\{d\}_k^T [A] \{d\}_k} $$
   2. Update the solution vector:
      $$ \{x\}_{k+1} = \{x\}_k + \alpha_k \{d\}_k $$
   3. Update the residual vector:
      $$ \{r\}_{k+1} = \{r\}_k - \alpha_k [A] \{d\}_k $$
   4. Convergence Check: If the magnitude of the residual \( ||\{r\}_{k+1}|| \) is less than a specified tolerance \( \epsilon \), exit the loop.
   5. Compute the improvement factor \( \beta_k \) to find the new conjugate direction:
      $$ \beta_k = \frac{\{r\}_{k+1}^T \{r\}_{k+1}}{\{r\}_k^T \{r\}_k} $$
   6. Update the search direction vector (ensuring it is \(A\)-conjugate to previous directions):
      $$ \{d\}_{k+1} = \{r\}_{k+1} + \beta_k \{d\}_k $$

Why use Conjugate Gradient over Gaussian Methods?

Gaussian elimination (a direct method) requires \(O(n^3)\) mathematical operations and structurally alters the matrix \([A]\) via row operations, a phenomenon known as "fill-in." This fill-in completely destroys the sparse nature of the matrix, requiring massive amounts of computer memory.

The Conjugate Gradient Method is vastly superior for massive, sparse systems because:

• No Fill-in: It relies solely on matrix-vector multiplications (\([A]\{d\}\)). The matrix \([A]\) is never altered, completely preserving its sparsity and minimizing memory requirements.
• Early Stopping: Being iterative, the algorithm can be halted as soon as the error (residual) falls below a practical engineering tolerance, providing a fast approximate solution without completing all \(n\) steps.

Limitations of the Conjugate Gradient Method

• Matrix Dependency: The matrix \([A]\) must strictly be symmetric and positive-definite. It fails on non-symmetric matrices.
• Conditioning: The rate of convergence is highly sensitive to the condition number of the matrix. Ill-conditioned matrices suffer from round-off errors and converge very slowly unless a "Preconditioned" Conjugate Gradient (PCG) method is utilized.

2.5 FOURIER TRANSFORM

Time Domain vs. Frequency Domain

In civil engineering dynamics (like earthquake engineering or wind vibrations), structural responses are recorded via sensors.

• Time Domain: Represents how the signal's amplitude changes over time. An earthquake accelerogram is a time-domain signal. Parameters identified: Peak Ground Acceleration (PGA), duration of strong motion, and arrival time of waves.
• Frequency Domain: Represents the signal broken down into its constituent sine/cosine waves at different frequencies. Parameters identified: The predominant frequency of the earthquake, the natural resonant frequencies of the building, and the amplitude (energy) associated with each harmonic.

Discrete Fourier Transform (DFT) and Fast Fourier Transform (FFT)

The Discrete Fourier Transform (DFT) is the mathematical algorithm used to convert a finite sequence of equally-spaced discrete time-domain samples into a same-length sequence of frequency-domain complex numbers. The formula is:

Computing the DFT directly takes \(O(N^2)\) operations, which is prohibitively slow for long sensor data files.

The Fast Fourier Transform (FFT), commonly implementing the Cooley-Tukey algorithm, is a highly optimized computational technique that computes the exact same result as the DFT but factorizes the matrix mathematically to reduce the computation time to \(O(N \log N)\). This makes real-time digital signal processing possible.

MATLAB Program for Fast Fourier Transform

Below is a standard algorithm implemented via a MATLAB script to perform an FFT on a given time-domain signal and plot the frequency spectrum.