Problem 2 - A First Course in the Numerical Analysis of Differential Equations

Question

Apply the plain conjugate gradient method to the linear system $$ \left[\begin{array}{lll} 1 & 0 & 0 \\ 0 & 2 & 0 \\ 0 & 0 & 3 \end{array}\right] \boldsymbol{x}=\left[\begin{array}{l} 1 \\ 1 \\ 1 \end{array}\right] $$ starting as usual with \(\boldsymbol{x}^{[0]}=\boldsymbol{0}\). Verify that the residuals \(\boldsymbol{r}^{[0]}, \boldsymbol{r}^{[1]}\) and \(\boldsymbol{r}^{[2]}\) are mutually orthogonal, that the search directions \(\boldsymbol{d}^{[0]}, \boldsymbol{d}^{[1]}\) and \(\boldsymbol{d}^{[2]}\) are mutually conjugate and that \(\boldsymbol{x}^{[3]}\) satisfies the linear system.

Step-by-Step Solution

Answer

The solution is \(\boldsymbol{x}^{[3]} = \begin{bmatrix} 1 \\ 1 \\ 1 \end{bmatrix}\). Residuals are mutually orthogonal; search directions are mutually conjugate.

Step 1: Define the Problem

We have a linear system of equations given by a diagonal matrix \(A = \begin{bmatrix} 1 & 0 & 0 \ 0 & 2 & 0 \ 0 & 0 & 3 \end{bmatrix}\) and a vector \(\boldsymbol{b} = \begin{bmatrix} 1 \ 1 \ 1 \end{bmatrix}\). We need to solve it using the conjugate gradient method starting with an initial guess \(\boldsymbol{x}^{[0]} = \boldsymbol{0}\).

Step 2: Initialize Variables

The initial guess for \(\boldsymbol{x}\) is \(\boldsymbol{x}^{[0]} = \begin{bmatrix} 0 \ 0 \ 0 \end{bmatrix}\). Compute the initial residual \(\boldsymbol{r}^{[0]} = \boldsymbol{b} - A\boldsymbol{x}^{[0]} = \boldsymbol{b} = \begin{bmatrix} 1 \ 1 \ 1 \end{bmatrix}\). Set the initial search direction \(\boldsymbol{d}^{[0]} = \boldsymbol{r}^{[0]} = \begin{bmatrix} 1 \ 1 \ 1 \end{bmatrix}\).

Step 3: Iteration 1

Calculate \(\alpha^{[0]} = \frac{(\boldsymbol{r}^{[0]})^T \boldsymbol{r}^{[0]}}{(\boldsymbol{d}^{[0]})^T A \boldsymbol{d}^{[0]}} = \frac{3}{6} = \frac{1}{2}\). Update estimate \(\boldsymbol{x}^{[1]} = \boldsymbol{x}^{[0]} + \alpha^{[0]} \boldsymbol{d}^{[0]} = \begin{bmatrix} \frac{1}{2} \ \frac{1}{2} \ \frac{1}{2} \end{bmatrix}\). Calculate \(\boldsymbol{r}^{[1]} = \boldsymbol{r}^{[0]} - \alpha^{[0]} A\boldsymbol{d}^{[0]} = \begin{bmatrix} 0 \ \frac{1}{2} \ 0 \end{bmatrix}\).

Step 4: Update Search Direction 1

Calculate \(\beta^{[0]} = \frac{(\boldsymbol{r}^{[1]})^T \boldsymbol{r}^{[1]}}{(\boldsymbol{r}^{[0]})^T \boldsymbol{r}^{[0]}} = \frac{1}{3}\). Update search direction \(\boldsymbol{d}^{[1]} = \boldsymbol{r}^{[1]} + \beta^{[0]} \boldsymbol{d}^{[0]} = \begin{bmatrix} \frac{1}{3} \ \frac{5}{6} \ \frac{1}{3} \end{bmatrix}\).

Step 5: Iteration 2

Calculate \(\alpha^{[1]} = \frac{(\boldsymbol{r}^{[1]})^T \boldsymbol{r}^{[1]}}{(\boldsymbol{d}^{[1]})^T A \boldsymbol{d}^{[1]}} = \frac{1}{5}\). Update \(\boldsymbol{x}^{[2]} = \boldsymbol{x}^{[1]} + \alpha^{[1]} \boldsymbol{d}^{[1]} = \begin{bmatrix} \frac{3}{5} \ \frac{3}{5} \ \frac{3}{5} \end{bmatrix}\). Calculate \(\boldsymbol{r}^{[2]} = \boldsymbol{r}^{[1]} - \alpha^{[1]} A\boldsymbol{d}^{[1]} = \begin{bmatrix} 0 \ 0 \ \frac{2}{5} \end{bmatrix}\).

Step 6: Verify Orthogonality

Check orthogonality of residuals: \(\boldsymbol{r}^{[0]} \cdot \boldsymbol{r}^{[1]} = 0\), \(\boldsymbol{r}^{[1]} \cdot \boldsymbol{r}^{[2]} = 0\), \(\boldsymbol{r}^{[0]} \cdot \boldsymbol{r}^{[2]} = 0\). All products are zero, confirming orthogonality.

Step 7: Final Iteration

Given mutually orthogonal residuals, calculate \(\alpha^{[2]} = 1\) to resolve entirely for \(\boldsymbol{x}^{[3]}\). Hence, \(\boldsymbol{x}^{[3]} = \begin{bmatrix} 1 \ 1 \ 1 \end{bmatrix}\), which satisfies the linear system \(A\boldsymbol{x} = \boldsymbol{b}\).

Step 8: Confirm Conjugate Directions

Verify that directions \(\boldsymbol{d}^{[0]}, \boldsymbol{d}^{[1]},\) and \(\boldsymbol{d}^{[2]}\) are conjugate by ensuring \((\boldsymbol{d}^{[i]})^T A \boldsymbol{d}^{[j]} = 0\) for \(i eq j\). This confirms that each search direction is indeed conjugate with respect to the matrix \(A\).

Key Concepts

Numerical Linear Algebra

Numerical Linear Algebra is the study of implementing matrix operations and algorithms on computers. In this exercise, we work with linear systems and the conjugate gradient method, a valuable tool for solving large symmetric positive-definite systems. The conjugate gradient method is iterative, which means it repeatedly improves an estimate to converge to a solution. By starting with an initial guess, here \(\boldsymbol{x}^{[0]} = \boldsymbol{0}\), the method updates this estimate step by step.
Some key computations involved are the residuals \(\boldsymbol{r}^{[i]}\), which reflect how far our current solution is from solving the linear system. In each iteration, the residual also helps determine how we update our guess. A core reason for the method's efficacy is its relation to minimizing a quadratic form given by the matrix. This is distinct from directly inverting the matrix, making the method efficient and effective even for large matrices.

• Advantages: Handles large systems effectively and doesn't require explicit matrix inversion.
• Applications: Useful in engineering, physics simulations, and large data processing tasks.

Orthogonality of Residuals

Orthogonality is a property where vectors stand at a right angle to each other, leading to a scalar dot product of zero. In the conjugate gradient method, ensuring that the residuals \(\boldsymbol{r}^{[i]}\) are orthogonal is crucial for convergence. When we say residuals are orthogonal, it means the remaining error in our solution from one iteration is independent of the errors in prior iterations.
This orthogonality guarantees that each step in the conjugate gradient method effectively reduces the error. In our exercise, we see that the calculated dot products \(\boldsymbol{r}^{[0]} \cdot \boldsymbol{r}^{[1]} = 0\), \(\boldsymbol{r}^{[1]} \cdot \boldsymbol{r}^{[2]} = 0\), and \(\boldsymbol{r}^{[0]} \cdot \boldsymbol{r}^{[2]} = 0\) confirm this crucial aspect of orthogonality. Successfully achieving orthogonal residuals indicates our solution estimates are improving without redundant adjustments.

• Purpose: Helps ensure distinct updates towards the solution.
• Verification: Checking dot products of consecutive residuals to ensure they are zero.

Mutual Conjugacy of Directions

The concept of mutual conjugacy pertains to the search directions, or how the estimates are updated during each step in the conjugate gradient method. Search directions \(\boldsymbol{d}^{[i]}\) are computed such that they are conjugate with respect to the matrix \(A\). Two vectors \(\boldsymbol{d}^{[i]}\) and \(\boldsymbol{d}^{[j]}\) are conjugate if their product with the matrix gives zero: \(\boldsymbol{d}^{[i]^T} A \boldsymbol{d}^{[j]} = 0\).
Establishing conjugacy ensures that each new direction chosen explores a fresh dimension of the matrix’s behavior, allowing for rapid convergence. In our step-by-step example, we ensure this mutual conjugacy by computing new directions using a combination of current residuals and previous directions. This strategic approach helps diminish overlap in the new directions, similar to how orthogonality does for residuals.

• Benefit: Ensures more efficient exploration of solution space.
• Check: Calculating \(\boldsymbol{d}^{[i]^T} A \boldsymbol{d}^{[j]} = 0\) helps verify mutual conjugacy.