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

Question

This question is concerned with the relevance of the linear stability domain to the numerical solution of inhomogeneous linear systems. a Let \(\Lambda\) be a nonsingular matrix. Prove that the solution of \(\boldsymbol{y}^{\prime}=\Lambda \boldsymbol{y}+\boldsymbol{a}\) \(\boldsymbol{y}\left(t_{0}\right)=\boldsymbol{y}_{0}\), is $$ y(t)=\mathrm{e}^{\left(t-t_{0}\right) \Lambda} \boldsymbol{y}_{0}+\Lambda^{-1}\left[\mathrm{e}^{\left(t-t_{0}\right) \Lambda}-I\right] \boldsymbol{a}, \quad t \geq t_{0} $$ Thus, deduce that if \(\Lambda\) has a full set of eigenvectors and all its the eigenvalues reside in \(\mathrm{C}^{-}\)then \(\lim _{t \rightarrow \infty} \boldsymbol{y}(t)=-\Lambda^{-1} \boldsymbol{a}\). b Assuming for simplicity's sake that the underlying equation is scalar, i.e. \(y^{\prime}=\lambda y+a, y\left(t_{0}\right)=y_{0}\), prove that a single step of the Runge-Kutta method (3.9) results in $$ y_{n+1}=r(h \lambda) y_{n}+q(h \lambda), \quad n=0,1, \ldots $$ where \(r\) is given by (4.13) and $$ q(z):=h a \boldsymbol{b}^{T}(I-z A)^{-1} \boldsymbol{c} \in \mathbb{P}_{(\nu-1) / \nu}, \quad z \in \mathrm{C} $$ c Deduce, by induction or otherwise, that $$ y_{n}=[r(h \lambda)]^{n} y_{0}+\left\\{\frac{[r(h \lambda)]^{n}-1}{r(h \lambda)-1}\right\\} q(h \lambda), \quad n=0,1, \ldots $$ d Assuming that \(h \lambda \in \mathcal{D}\), prove that \(\lim _{n \rightarrow \infty} y_{n}\) exists and is bounded.

Step-by-Step Solution

Answer

The solution exists and is bounded as \(n \to \infty\) if \(|r(h \lambda)| < 1\).

Step 1: Understand the solution to the differential equation

The given differential equation is \(\boldsymbol{y}^{\prime}=\Lambda \boldsymbol{y}+\boldsymbol{a}\) with initial condition \(\boldsymbol{y}(t_0)=\boldsymbol{y}_0\). Its solution can be written using the matrix exponential. The general solution is given as \(y(t)=\mathrm{e}^{(t-t_0) \Lambda}\boldsymbol{y}_0+\Lambda^{-1}[\mathrm{e}^{(t-t_0) \Lambda}-I] \boldsymbol{a}\). This is derived by integrating the homogeneous and particular solutions.

Step 2: Analyze the asymptotic behavior

To determine the behavior as \(t \to \infty\), note that if \(\Lambda\) has a full set of eigenvectors and eigenvalues in \(\mathbb{C}^{-}\), \(\mathrm{e}^{(t-t_0)\Lambda} \to 0\) as \(t \to \infty\). Hence, \(\lim_{t \to \infty} \boldsymbol{y}(t) = -\Lambda^{-1} \boldsymbol{a}\), as the exponential term vanishes.

Step 3: Express Runge-Kutta for scalar case

For the scalar case, the Runge-Kutta method is applied to \(y' = \lambda y + a\) with \(y(t_0) = y_0\). The result after a single step is expressed as \(y_{n+1} = r(h \lambda) y_n + q(h \lambda)\), where \(r\) and \(q\) come from the method's update rules.

Step 4: Prove the iterative process result

By induction or recursive substitution, \(y_n = [r(h \lambda)]^n y_0 + \{\frac{[r(h \lambda)]^n - 1}{r(h \lambda) - 1}\} q(h \lambda)\) can be established. The base case is \(n=0\), and the inductive step uses the form of \(y_{n+1}\) in Step 3.

Step 5: Deduce the limit behavior of \\(y_n\\)

Assuming \(h \lambda\) is in the linear stability domain \(\mathcal{D}\), where \(|r(h \lambda)| < 1\), it follows that \([r(h \lambda)]^n \to 0\) as \(n \to \infty\). Thus, \(\lim_{n \to \infty} y_n = \lim_{n \to \infty} \left\{ \frac{[r(h \lambda)]^n - 1}{r(h \lambda) - 1} \right\} q(h \lambda) = \frac{-1}{r(h \lambda) - 1} q(h \lambda)\), which is bounded.

Key Concepts

Linear Stability

Linear stability is an essential concept when solving differential equations numerically, ensuring that the numerical solution behaves similarly to the true solution. It is especially crucial in iterative methods, where errors can grow rapidly if the system is unstable. Linear stability refers to a property of a numerical method that determines whether small perturbations to the solution will grow or decay over time. If these perturbations decay, the system is considered stable; otherwise, it is unstable.

In numerical analysis, we often explore the stability region of a numerical method, specifically within the complex plane. A method is said to have linear stability if, for each eigenvalue of the matrix governing the system, the corresponding point in the complex plane lies within the stability region. This implies that the method can handle various components of the solution's behavior without the risk of numerical explosions. Understanding linear stability helps in choosing the right numerical method and step sizes for a given problem.

Runge-Kutta Method

The Runge-Kutta method is a widely used numerical technique for solving ordinary differential equations (ODEs). It enhances accuracy compared to simpler methods like the Euler method, without requiring complex calculations at each step. The Runge-Kutta method updates the solution using weighted averages of slope estimates at various points, typically involving multiple evaluations within a single time step.

In the context of the given exercise, the Runge-Kutta method is applied to solve a first-order differential equation of the form \(y' = \lambda y + a\). The resulting formula after a single step can be expressed as \(y_{n+1} = r(h \lambda) y_n + q(h \lambda)\). Here, \(r\) and \(q\) are polynomials derived from the method's coefficients, providing a systematic way to update the solution iteratively. This method ensures stability when correctly implemented, maintaining both efficiency and accuracy.

Matrix Exponential

The matrix exponential is a mathematical function that extends the concept of exponential functions to matrices. It is vital in solving systems of linear differential equations. For a matrix \( \Lambda \), the matrix exponential \( \mathrm{e}^{t \Lambda} \) gives the fundamental solution matrix to the homogeneous linear system \( \mathbf{y}' = \Lambda \mathbf{y} \).

In the exercise provided, the matrix exponential plays a crucial role in expressing the solution for the inhomogeneous system \( \mathbf{y}' = \Lambda \mathbf{y} + \mathbf{a} \). The formula \(y(t)=\mathrm{e}^{(t-t_{0}) \Lambda}\boldsymbol{y}_{0}+\Lambda^{-1}[\mathrm{e}^{(t-t_{0}) \Lambda}-I] \boldsymbol{a}\) demonstrates how the matrix exponential aids in computing both the transient and steady-state behaviors. It accounts for the initial condition and any forced component \( \mathbf{a} \). Understanding the matrix exponential's properties, such as its relationship to eigenvalues, is essential for working with linear differential systems.

Eigenvalues

Eigenvalues are critical to understanding the behavior of systems represented by matrices, especially in the context of stability and dynamics. An eigenvalue of a matrix is a scalar that provides insight into the matrix's directional scaling properties. When solving linear differential equations, eigenvalues are essential in analyzing the system's response over time.

In the scope of the exercise, the eigenvalues of matrix \( \Lambda \) help determine the asymptotic behavior of the system \( \mathbf{y}' = \Lambda \mathbf{y} + \mathbf{a} \). If the eigenvalues have negative real parts, the exponentials decay, leading the system to tend toward a steady-state. The limit \( \lim_{t \rightarrow \infty} \mathbf{y}(t) = - \Lambda^{-1} \mathbf{a} \) holds true under this condition, demonstrating that eigenvalues are a decisive factor in predicting long-term outcomes and ensuring stability in numerical simulations.