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

Question

We consider the one-step stationary method \(M \boldsymbol{x}^{[k+1]}=N \boldsymbol{x}^{[k]}+\boldsymbol{b}\), where \(M-N=A\) and the matrix \(M\) is nonsingular. a Prove that \(\boldsymbol{x}^{[k]}-\boldsymbol{x}=H^{k} \boldsymbol{e}^{[0]}\), where \(H=M^{-1} N\) is the iteration matrix, \(e^{[0]}=x^{[0]}-x\) and \(x\) is the exact solution of the linear system. b Given \(m \geq 1\), we form a new candidate solution \(y^{[m]}\) by the linear combination $$ \boldsymbol{y}^{[m]}=\sum_{k=0}^{m} \nu_{k} \boldsymbol{x}^{[k]}, \quad \text { where } \quad \sum_{k=0}^{m} \nu_{k}=1 $$ Prove that \(\boldsymbol{y}^{[m]}-\boldsymbol{x}=\sum_{k=0}^{m} \nu_{k} H^{k} \boldsymbol{e}^{[0]}\), and thus deduce that $$ \left\|\boldsymbol{y}^{[m]}-\boldsymbol{x}\right\|_{2} \leq\left\|p_{m}(H)\right\|_{2}\left\|\boldsymbol{e}^{[0]}\right\|_{2}, \quad \text { where } \quad p(z)=\sum_{k=0}^{m} \nu_{k} z^{k} $$ c Suppose that it is known that all the eigenvalues of \(H\) are real and reside in the interval \([\alpha, \beta]\) and that the matrix has a full set of eigenvectors. Prove that $$ \left\|\boldsymbol{y}^{[m]}-\boldsymbol{x}\right\|_{2} \leq\|V\|_{2}\left\|V^{-1}\right\|_{2}\left\|\boldsymbol{e}^{[0]}\right\|_{2} \max _{x \in[\alpha, \beta]}\left|p_{m}(x)\right| $$ where \(V\) is the matrix of the eigenvectors of \(H\). d We now use our freedom of choice of the parameters \(\nu_{0}, \ldots, \nu_{m}\), hence of the polynomial \(p_{m}\) such that \(p_{m}(1)=1\) (why this condition?), to minimize \(|p(x)|\) for \(x \in[\alpha, \beta]\). To this end prove that the Chebyshev polynomial \(T_{m}\) (see Exercise \(3.2\) ) satisfies the inequality \(\left|T_{m}(x)\right| \leq 1, x \in[-1,1]\). (Since \(T_{n}(1)=1\), this inequality cannot be improved by any other polynomial \(q\) such that \(\max _{x \in[-1,1]}|q(x)|=1\).) Deduce that the best choice of \(p_{m}\) is $$ p_{m}(x)=\frac{T_{m}(2(x-\alpha) /(\beta-\alpha)-1)}{T_{m}(2(1-\alpha) /(\beta-\alpha)-1)} . $$ e Show that this algorithm can be formulated in a Krylov subspace formalism. (This is the famed Chebyshev iterative method. Note, however, that naive implementation of this iterative procedure is problematic.)

Step-by-Step Solution

Answer

Use iterative steps: prove relations, form candidates, bound with eigenvalues, minimize with Chebyshev polynomials in Krylov subspace form.

Step 1: Show the relation for the one-step method

The relation given is \(M \boldsymbol{x}^{[k+1]}=N \boldsymbol{x}^{[k]} + \boldsymbol{b}\). Rearrange this to get \( \boldsymbol{x}^{[k+1]} = M^{-1}N \boldsymbol{x}^{[k]} + M^{-1}\boldsymbol{b}\). For the exact solution \( \boldsymbol{x} \), it must satisfy \( A \boldsymbol{x} = \boldsymbol{b} \). Thus, \( \boldsymbol{x} = M^{-1}N\boldsymbol{x} + M^{-1}\boldsymbol{b} \). By subtracting these, we get \(\boldsymbol{x}^{[k+1]} - \boldsymbol{x} = M^{-1}N\boldsymbol{x}^{[k]} - M^{-1}N\boldsymbol{x} = M^{-1}N(\boldsymbol{x}^{[k]} - \boldsymbol{x})\). Therefore, our iteration matrix \(H\) is \( H = M^{-1}N \) and we prove \( \boldsymbol{x}^{[k]} - \boldsymbol{x} = H^{k}\boldsymbol{e}^{[0]} \) recursively.

Step 2: Formulate a new candidate with coefficients

Given a new candidate solution \(\boldsymbol{y}^{[m]} = \sum_{k=0}^{m} u_k \boldsymbol{x}^{[k]} \), using step 1 \(\boldsymbol{x}^{[k]} - \boldsymbol{x} = H^{k} \boldsymbol{e}^{[0]} \), substitute into the expression for \( \boldsymbol{y}^{[m]} \) to get: \(\boldsymbol{y}^{[m]} - \boldsymbol{x} = \sum_{k=0}^{m} u_k (\boldsymbol{x}^{[k]} - \boldsymbol{x}) = \sum_{k=0}^{m} u_k H^k \boldsymbol{e}^{[0]}\). Since \(\sum_{k=0}^{m} u_k = 1\), we obtain \(\left\|\boldsymbol{y}^{[m]}-\boldsymbol{x}\right\|_{2} \leq \left\|p_{m}(H)\right\|_{2} \left\|\boldsymbol{e}^{[0]}\right\|_{2}\), where \(p_m(z) = \sum_{k=0}^{m} u_k z^k\).

Step 3: Analyze bounding with real eigenvalues

Given eigenvalues of \(H\) are real within \([\alpha, \beta]\), use \(\boldsymbol{y}^{[m]} - \boldsymbol{x} = \sum_{k=0}^{m} u_k H^k \boldsymbol{e}^{[0]}\). Rewrite using a basis of eigenvectors \(V\) such that \(H = V\Lambda V^{-1}\), resulting in \(p_{m}(H) = V p_{m}(\Lambda) V^{-1}\). Thus, \(\left\| \boldsymbol{y}^{[m]} - \boldsymbol{x} \right\|_{2} \leq \|V\|_{2}\left\|V^{-1}\right\|_{2} \left\|\boldsymbol{e}^{[0]}\right\|_{2} \max_{x \in [\alpha, \beta]} |p_{m}(x)|\).

Step 4: Use Chebyshev polynomials for minimization

The Chebyshev polynomial \(T_m\) satisfies \(\left|T_m(x)\right| \leq 1\) for \(x \in [-1, 1]\). To minimize \(|p(x)|\) for \(x \in [\alpha, \beta]\) under the condition \(p_m(1) = 1\), construct \(p_m(x) = \frac{T_m(2(x-\alpha)/ (\beta-\alpha) - 1)}{T_m(2(1-\alpha)/(\beta-\alpha) - 1)}\). Use translation and scaling based on \([\alpha, \beta]\) to fit inside \([-1,1]\).

Step 5: Connect to Krylov subspace

The algorithm reduces the error over iterations and optimizes it through a sequence of polynomial approximations using Chebyshev polynomials, representing a Krylov subspace method. The basis vectors of this subspace correspond to the powers of \(H\) acting on \(\boldsymbol{e}^{[0]}\), suitable for iterative methods as in Krylov subspaces, allowing the representation as sums of powers influenced by the initial vector.

Key Concepts

Numerical Analysis

Numerical analysis is a branch of mathematics that focuses on developing algorithms to solve mathematical problems numerically, as opposed to symbolically. This is incredibly useful when problems do not have exact analytical solutions or when calculations must be performed with a high degree of precision. Numerical methods are applicable across numerous fields such as physics, engineering, and finance, where approximations are often necessary due to the complexity of models or the sheer volume of data involved.

Key aspects of numerical analysis include:

• Stability: Ensuring that small changes in input or rounding errors don't lead to significant variations in output.
• Convergence: Ensuring that an iterative method approaches the exact solution as the number of iterations increases.
• Efficiency: Assessing the computational resources required by the algorithm to reach a solution.

In the context of the Chebyshev iterative method, numerical analysis helps us understand and control how approximate solutions to linear systems are derived step by step. It ensures that the iterative processes yield results that are numerically stable and converge efficiently to the true solution. Practical implementations often prioritize such algorithms to balance precision, speed, and resource use.

Iterative Methods

Iterative methods are techniques used in numerical analysis to solve mathematical problems by finding successive approximations to the solutions. Unlike direct methods, which attempt to solve problems in a finite number of steps (like solving a system of linear equations using Gaussian elimination), iterative methods start with an initial guess and refine it through repeated applications of an algorithm.

Some advantages of iterative methods include:

• Flexibility: They can handle large, sparse systems effectively, which are common in simulations and real-world problems.
• Adaptability: These methods are useful when the system is too large for direct methods, allowing for parallel computations.
• Resource efficiency: By improving estimates incrementally, they often require less memory.

The Chebyshev iterative method fits into this framework by iteratively refining a solution approximation through polynomial transformations. It leverages properties of Chebyshev polynomials to minimize errors efficiently. This method is particularly powerful when dealing with linear systems where the iterative matrix's spectral properties can be exploited to speed up convergence.

Eigenvalues and Eigenvectors

Understanding eigenvalues and eigenvectors is crucial when working with many numerical methods, including the Chebyshev iterative method. For a given square matrix, an eigenvalue is a scalar that indicates how much a corresponding eigenvector is scaled during a transformation. If you multiply a matrix by a vector and get a resultant vector that is just a scaled version of the original, that vector is an eigenvector, and the scaling factor is the eigenvalue.

Some important features include:

• Eigenvalues as Scalars: They define the characteristic roots of the matrix and indicate stability properties.
• Eigenvectors as Bases: These can form a basis that simplifies complex matrix operations, aiding in the analysis and solutions of systems.
• Spectral Theorem: If a matrix can be diagonalized, its behavior can be understood through its eigenvalues and eigenvectors, simplifying many problems.

In the Chebyshev iterative method, the eigenvalues of the iteration matrix determine the convergence behavior. The method effectively minimizes errors by adjusting polynomial interpolations based on these eigenvalues, ensuring efficient convergence across real intervals.

Stationary Methods

Stationary methods are a subset of iterative methods used for solving linear algebraic systems where the algorithm is fixed at each iteration. These methods do not adjust their parameters from one iteration to the next. A common example is the Jacobi or Gauss-Seidel methods.

The properties of stationary methods include:

• Simplicity: They are defined by a fixed procedure that doesn't change over iterations.
• Predictability: Due to their non-varying nature, their performance can be easy to analyze theoretically.
• Foundation for Methods: They are often used as foundational building blocks for more advanced iterative techniques.

In the exercise at hand, the Chebyshev iterative method can be considered a stationary method as it employs a fixed matrix, derived from the original problem, across iterations. However, it cleverly enriches this fixed approach by incorporating Chebyshev polynomials to optimize convergence for a specified interval, enhancing the overall efficiency and robustness of solution processes.