Consistency, Stability, and Convergence

(20 pts) For the ODE \(u'(t) = f(t,u)\), consider the time-integration scheme \[ u_{n+1} = u_n + \frac {k}{4}f(t_n,u_n) + \frac {3k}{4}f\left (t_n+\frac {2}{3}k,\ u_n+\frac {2}{3}kf(t_n,u_n)\right ). \]

1. (5 pts) Identify the Butcher tableau for this scheme.
2. (8 pts) Compute the order of consistency of this scheme.
3. (7 pts) Compute the stability/amplification factor of the scheme.

(20 pts)

1. (10 pts) Use von Neumann stability analysis to determine a stability condition for, \[ D^+u_j^n = D_0u_j^n, \] where \(D^+u_j^n = (u_j^{n+1}-u_j^n)/k\), \(D_0u_j^n = (u_{j+1}^n-u_{j-1}^n)/(2h)\), and \(D_\pm u_j^n = \pm (u_{j\pm 1}^n-u_j^n)/h\).
2. (10 pts) Use von Neumann stability analysis to determine a stability condition for \[ D^+u_j^n = D_+D_-u_j^n. \]

(20 pts) Consider the PDE, \[ u_t + au_x = 0,\qquad u(x,0)=u_0(x),\qquad a > 0, \] with periodic boundary conditions over the spatial domain \(x \in [0, 1)\). Use the standard notation, \(u_j^n \approx u(x_j,t^n)\) with \(x_j = jh\) and \(t^n = nk\) with \(h = \frac {1}{M} > 0\) for some \(M \in \mathbb {N}\) and \(k > 0\). Consider the following semi-Lagrangian scheme: \[ u_j^{n+1}=p_j(y),\qquad y := X(t^n;x_j,t^{n+1}), \] where \(X(t; x, t_0)\) is the time-\(t\) spatial location of the PDE’s characteristic curve passing through spatial location \(x\) at time \(t_0\), and \(p_j(\cdot )\) is the linear interpolant formed by the data \((x_{j-1},u_{j-1}^n)\) and \((x_j,u_j^n)\).

1. (10 pts) Write this scheme as an explicit function of the \(u_j^n\) values.
2. (10 pts) Use von Neumann stability analysis to determine a stability condition for this scheme.

(20 pts)

1. (10 pts) Compute the stability/amplification factor for forward Euler to solve \(u' = f(t,u)\), and plot the region of the stability.
2. (10 pts) Compute the stability/amplification factor for the trapezoidal rule/Crank- Nicolson scheme and use this to show that the scheme is A-stable.

(20 pts)

1. (10 pts) Use von Neumann stability analysis to determine a stability condition for, \[ D^+u_j^n = D_0u_j^n, \] where \(D^+u_j^n = (u_j^{n+1}-u_j^n)/k\), \(D_0u_j^n = (u_{j+1}^n-u_{j-1}^n)/(2h)\), and \(D_\pm u_j^n = \pm (u_{j\pm 1}^n-u_j^n)/h\).
2. (10 pts) Use von Neumann stability analysis to determine a stability condition for \[ D^+u_j^n = D_+D_-u_j^n. \]

(Linear Multistep Methods).

1. Define linear multistep method (give formula). Give definition of the region of absolute stability.
2. Show that the region of absolute stability for the trapezoidal method is the set of all complex \(h\lambda \) with \(\operatorname {Real}(\lambda ) < 0\).

( Heat Equation and Stability of the Scheme). Consider the implicit in time, Backward Euler method for the solution of the heat equation: \[ \frac {u_m^{n+1}-u_m^n}{\Delta t} -10\frac {u_{m+1}^{n+1}-2u_m^{n+1}+u_{m-1}^{n+1}}{h^2} = f_m^{n+1}, \] \[ u_m^0 = g_m,\qquad m = 0,\pm 1,\pm 2,\ldots ,\qquad n = 0,1,\ldots ,[T/\Delta t]-1, \]

and investigate the stability of the scheme using the von Neumann analysis.

Consider the boundary value problem \[ u''(x) = f(x),\quad \text {for }0 < x < 1, \] \[ u(0) = 0,\quad u(1) = 0, \]

where \(f\) is a smooth function on \([0, 1]\), and the finite difference scheme \[ \frac {U_{j-1}-2U_j+U_{j+1}}{h^2} = f(x_j), \qquad j = 1,2,\ldots ,m, \] for finding approximations \(U_j\) to \(u(x_j)\). Here, \(x_j = jh\) and \((m+1)h = 1\).

Analyze the accuracy, stability, and convergence properties of the scheme.

For the initial value problem \(u'(t) = f(u(t),t)\), \(u(0) = \eta \), where \(f(u,t)\) is continuous with respect to \(t\) and Lipschitz continuous with respect to \(u\), consider the scheme \[ \frac {U^{n+2}-\frac {4}{3}U^{n+1}+\frac {1}{3}U^n}{k} = \frac {2}{3}f(U^{n+2},t_{n+2}), \] where \(t_n = nk\), and \(U^n\) is supposed to approximate \(u(t_n)\).

1. Analyze the consistency, zero-stability, and convergence of this scheme.

2. If you apply this scheme to the initial value problem

   \[ u'(t) = -10^{12}(u(t)-\cos (t))-\sin (t),\qquad u(0) = 2, \]

   what issues should you consider in choosing the timestep \(k\)? Are there time intervals in which \(k\) must be very small to get a reasonable solution and others in which it does not? Explain your answers.

1. Consider the initial value problem for the constant-coefficient diffusion equation (with \(\beta > 0\)) \[ v_t = \beta v_{xx},\quad t > 0 \]

   with initial data \(v(x,0) = f(x)\). A scheme for this problem is: \[ \frac {u_j^{n+1}-u_j^n}{k} = \frac {\beta }{h^2}\left \{u_{j-1}^{n+1}-2u_j^{n+1}+u_{j+1}^{n+1}\right \}. \]

   Analyze the 2-norm stability of this scheme. For which values of \(k > 0\) and \(h > 0\) is the scheme stable? (Note that there are no boundary conditions here.)

2. Consider the variable coefficient diffusion equation \[ v_t = (\beta (x)v_x)_x,\quad 0 < x < 1,\quad t > 0 \]

   with Dirichlet boundary conditions \[ v(0,t) = 0,\qquad v(1,t) = 0 \]

   and initial data \(v(x,0) = f(x)\). Assume that \(\beta (x) \ge \beta _0 > 0\), and that \(\beta (x)\) is smooth. Let \(\beta _{j+1/2} = \beta (x_{j+1/2})\). A scheme for this problem is: \[ \frac {u_j^{n+1}-u_j^n}{k} = \frac {1}{h^2}\left \{\beta _{j-1/2}u_{j-1}^{n+1} -(\beta _{j-1/2}+\beta _{j+1/2})u_j^{n+1} +\beta _{j+1/2}u_{j+1}^{n+1}\right \}. \]

   Analyze the 2-norm stability of this scheme for solving this initial boundary value problem. DO NOT NEGLECT THE FACT THAT THERE ARE BOUNDARY CONDITIONS!

For the initial value problem \(u'(t) = f(u(t),t)\), \(u(0) = \eta \), where \(f(u,t)\) is continuous with respect to \(t\) and Lipschitz continuous with respect to \(u\), consider the scheme \[ \frac {U^{n+2}-\frac {4}{3}U^{n+1}+\frac {1}{3}U^n}{k} = \frac {2}{3}f(U^{n+2},t_{n+2}), \] where \(t_n = nk\), and \(U^n\) is supposed to approximate \(u(t_n)\).

1. Analyze the consistency, zero-stability, and convergence of this scheme.

2. If you apply this scheme to the initial value problem

   \[ u'(t) = -10^{12}(u(t)-\cos (t))-\sin (t),\qquad u(0) = 2, \]

   what issues should you consider in choosing the timestep \(k\)? Are there time intervals in which \(k\) must be very small to get a reasonable solution and others in which it does not? Explain your answers.

 Heat Equation and Stability of the Scheme

Consider the implicit in time, Backward Euler method for the solution of the heat equation: \[ \frac {u_m^{n+1}-u_m^n}{\Delta t} - 4\frac {u_{m+1}^{n+1}-2u_m^{n+1}+u_{m-1}^{n+1}}{h^2} = f_m^{n+1}, \] \[ u_m^0=g_m,\qquad m=0,\pm 1,\pm 2,\ldots ,\qquad n=0,1,\ldots ,[T/\Delta t]-1, \] and investigate the stability of the scheme using the von Neumann analysis.

Region of Absolute Stability and Example of Linear One-Step Method

1. Define a region of absolute stability for the linear multistep method.

2. Find the region of absolute stability for the Backward Euler Method.

The following function \(f(x)\) [\(x\) is a real variable] is expanded in a series in powers of \((x+3)\). Find the radius of convergence of the series in each case. \[ \begin {array}{rcl} \text {(a)}\quad f(x)&=&e^{x^2},\\[1.4ex] \text {(b)}\quad f(x)&=&\left (\dfrac {\sin x-3}{\sin x-2}\right )^2,\\[1.8ex] \text {(c)}\quad f(x)&=&\dfrac {\sin (x+1)}{(x+1)(x+3)}. \end {array} \]

Heat Equation Stability

Consider the variable coefficient diffusion equation \[ v_t=(\beta (x)v_x)_x,\qquad 0<x<1,\quad t>0 \] with Dirichlet boundary conditions \[ v(0,t)=0,\qquad v(1,t)=0 \] and initial data \(v(x,0)=f(x)\). Assume that \(\beta (x)\ge \beta _0>0\), and that \(\beta (x)\) is smooth. Let \(\beta _{j+1/2}=\beta (x_{j+1/2})\). A scheme for this problem is: \[ \frac {u_j^{n+1}-u_j^n}{k} = \frac {1}{h^2} \left \{ \beta _{j-1/2}u_{j-1}^{n+1} -(\beta _{j-1/2}+\beta _{j+1/2})u_j^{n+1} +\beta _{j+1/2}u_{j+1}^{n+1} \right \}. \] Analyze the 2-norm stability of this scheme for solving this initial boundary value problem.

DO NOT NEGLECT THE FACT THAT THE PROBLEM HAS VARIABLE COEFFICIENTS AND THAT THERE ARE BOUNDARY CONDITIONS AT 0 AND 1!

Function \(f(x)\) (see below) is expanded in powers of \(x-x_0\), i.e. in the Taylor or Laurent series with center \(x_0\) (\(x\) is a real variable, \(x_0\) is a real number). Find the radius of convergence of this series. \[ \begin {array}{rcll} \text {(a)}\quad f(x)&=&\dfrac {\sin x^2}{(x-3)^2(x^2+1)},& x_0=3,\\[1.4ex] \text {(b)}\quad f(x)&=&\dfrac {\sin x}{x},& x_0=5,\\[1.4ex] \text {(c)}\quad f(x)&=&\dfrac {\sin x}{\sin x+3},& x_0=8. \end {array} \]

Trapezoid rule

For the ODE \(y'=f(t,y)\), the Trapezoid rule is \[ \frac {y_{n+1}-y_n}{\Delta t} = \frac {1}{2}\left (f(t_n,y_n)+f(t_{n+1},y_{n+1})\right ). \] Compute the region of absolute stability for the Trapezoid rule. Is this method \(A\)-stable?

Unconditionally stable scheme

For the one-dimensional heat equation, \[ u_t=u_{xx}, \] consider the Crank-Nicolson, central difference scheme \[ \frac {u_j^{n+1}-u_j^n}{\Delta t} = \frac {1}{2} \left ( \frac {u_{j+1}^{n+1}-2u_j^{n+1}+u_{j-1}^{n+1}}{h^2} + \frac {u_{j+1}^{n}-2u_j^{n}+u_{j-1}^{n}}{h^2} \right ) \] where \(\Delta t\) and \(h\) are temporal and spatial stepsizes, respectively, and where \(u_j^n\) is the numerical approximation to \(u(x=jh,t=n\Delta t)\), with \(x_0=0\) and the initial time \(t=0\). Use von Neumann stability analysis to show that this scheme is unconditionally stable.

Central scheme

Consider the one-dimensional advection equation, \[ u_t=cu_x, \] where \(c\in \mathbb {R}\) is the wavespeed. Use von Neumann stability analysis to determine the values of \(c\) for which the following the central difference, forward Euler scheme is stable: \[ \frac {u_j^{n+1}-u_j^n}{\Delta t} = c\frac {u_{j+1}^{n}-u_{j-1}^{n}}{2h}. \]

Heat Equation and Stability of the Scheme

Consider the implicit in time, Backward Euler method for the solution of the heat equation: \[ \frac {u_m^{n+1}-u_m^n}{\Delta t} - \frac {u_{m+1}^{n+1}-2u_m^{n+1}+u_{m-1}^{n+1}}{h^2} = f_m^{n+1}, \] \[ u_m^0=g_m,\qquad m=0,\pm 1,\pm 2,\ldots ,\qquad n=0,1,\ldots ,[T/\Delta t]-1. \] and investigate the stability of the scheme using the von Neumann analysis.

Heat Equation Stability:

Consider the variable coefficient diffusion equation \[ v_t=(\beta (x)v_x)_x,\qquad 0<x<1,\quad t>0 \] with Dirichlet boundary conditions \[ v(0,t)=0,\qquad v(1,t)=0 \] and initial data \(v(x,0)=f(x)\). Assume that \(\beta (x)\ge \beta _0>0\), and that \(\beta (x)\) is smooth. Let \(\beta _{j+1/2}=\beta (x_{j+1/2})\). A scheme for this problem is: \[ \frac {u_j^{n+1}-u_j^n}{k} = \frac {1}{h^2} \left \{ \beta _{j-1/2}u_{j-1}^{n+1} -(\beta _{j-1/2}+\beta _{j+1/2})u_j^{n+1} +\beta _{j+1/2}u_{j+1}^{n+1} \right \}. \] Analyze the 2-norm stability of this scheme for solving this initial boundary value problem.

DO NOT NEGLECT THE FACT THAT THE PROBLEM HAS VARIABLE COEFFICIENTS AND THAT THERE ARE BOUNDARY CONDITIONS AT 0 AND 1!

Fact 1: A real symmetric \(n\times n\) matrix \(A\) can be diagonalized by an orthogonal similarity transformation, and \(A\)’s eigenvalues are real.

Fact 2: The \((N-1)\times (N-1)\) matrix \(M\) defined by \[ M= \begin {pmatrix} -2 & 1 & 0 & \cdots & 0\\ 1 & -2 & 1 & \ddots & \vdots \\ 0 & 1 & -2 & \ddots & 0\\ \vdots & \ddots & \ddots & \ddots & 1\\ 0 & \cdots & 0 & 1 & -2 \end {pmatrix} \] has eigenvalues \[ \mu _l=-4\sin ^2\left (\frac {\pi l}{2N}\right ), \qquad l=1,2,\ldots ,N-1. \]

Fact 3: The \((N+1)\times (N+1)\) matrix: \[ \begin {pmatrix} -1 & 1 & 0 & \cdots & 0\\ 1 & -2 & 1 & \ddots & \vdots \\ 0 & 1 & -2 & \ddots & 0\\ \vdots & \ddots & \ddots & \ddots & 1\\ 0 & \cdots & 0 & 1 & -1 \end {pmatrix} \] has eigenvalues \[ \mu _l=-4\sin ^2\left (\frac {\pi l}{2(N+1)}\right ), \qquad l=0,1,\ldots ,N. \]

Fact 4: For a real \(n\times n\) matrix \(A\), the Rayleigh quotient of a vector \(x\in \mathbb {R}^n\) is the scalar \[ r(x)=\frac {x^T Ax}{x^T x}. \] The gradient of \(r(x)\) is \[ \nabla r(x)=\frac {2}{x^T x}(Ax-r(x)x). \] If \(x\) is an eigenvector of \(A\) then \(r(x)\) is the corresponding eigenvalue and \(\nabla r(x)=0\).

The following function \(f(x)\) [\(x\) is a real variable] can be represented by a series (Taylor or Laurent) in powers of \(x\). Find the radius of convergence of the series in each case \[ \begin {array}{rcl} \text {(a)}\quad f(x)&=&\left (\dfrac {\sin x}{x}\right )^2,\\[1.4ex] \text {(b)}\quad f(x)&=&\dfrac {\sin x}{x^2},\\[1.4ex] \text {(c)}\quad f(x)&=&\left (\dfrac {x}{\sin x}\right )^2,\\[1.4ex] \text {(d)}\quad f(x)&=&\dfrac {x^2}{\sin x},\\[1.4ex] \text {(e)}\quad f(x)&=&\dfrac {\sin x}{x^2+1}. \end {array} \] [You do not need to find the series themselves.]

The following function \(f(x)\) [\(x\) is a real variable] can be represented by a series (Taylor or Laurent) in powers of \((x-7)\). Find the radius of convergence of the series in each case \[ \begin {array}{rcl} \text {(a)}\quad f(x)&=&e^{(x-7)^{10}},\\[0.6ex] \text {(b)}\quad f(x)&=&\left (\dfrac {\sin x-3}{\sin x-2}\right )^2,\\[1.4ex] \text {(c)}\quad f(x)&=&\dfrac {\sin x-2}{x^2-49}. \end {array} \] [You do not need to find the series themselves.]

Linear Multistep Methods

Give the definition of linear multistep method (give formula).

a) State necessary and sufficient condition for linear multistep method to be consistent.

b) State root condition for linear multistep method.

c) Construct an example of a consistent but not stable linear multistep method. Justify your answer.

Region of Absolute Stability and Example of Linear One-Step Method

1. Define a region of absolute stability for Linear Multistep Method.

2. Find the region of absolute stability for the Forward Euler Method.

Upwind Scheme

Consider the advection equation \[ u_t+5u_x=0,\qquad x_L<x<x_R,\qquad 0<t\le T, \] where \(u(x,0)=g(x)\), and \(u(x_L,t)=u_L(t)\) for \(t>0\).

a) Write the Upwind Scheme for this problem.

b) What is the local truncation error of the Upwind Scheme?

c) What is the stencil of the scheme? What is the CFL condition for this method?

d) Investigate the stability of the method using Von Neumann Stability Analysis.

Linear Multistep Methods

Define linear multistep method (give formula).

a) State necessary and sufficient condition for linear multistep method to be consistent.

b) State root condition for linear multistep method.

c) What can you say about consistency and stability properties of the method: \[ y_{n+1}-3y_n+2y_{n-1}=-h f(y_{n-1})? \]

Upwind Scheme

Consider the advection equation \[ u_t-u_x=0,\qquad x_L<x<x_R,\qquad 0<t\le T, \] where \(u(x,0)=g(x)\), and \(u(x_R,t)=u_R(t)\) for \(t>0\).

a) Write the Upwind Scheme for this problem.

b) What is the stencil of the scheme? What is the CFL condition for this method?

c) Investigate the stability of the method using Von Neumann Stability Analysis.

The following functions \(f(x)\) [\(x\) is a real variable] can be represented by series (Taylor or Laurent) in powers of \((x-2)\). Find the radius of convergence of the series in each case \[ \begin {array}{rcl} \text {(a)}\quad f(x)&=&e^{x^2},\\[0.6ex] \text {(b)}\quad f(x)&=&\dfrac {\sin x-3}{\sin x-2},\\[1.2ex] \text {(c)}\quad f(x)&=&\dfrac {e^x}{(x-2)(1+e^x)}. \end {array} \]

ODEs Consider the initial value problem for a system of ODEs \(y'=Ay\) where \(A\) has eigenvalues that are widely spread in the left half plane. Which of the following schemes would be the best choice for solving this problem? Justify your answer in terms of stability, accuracy, and efficiency. \[ \begin {aligned} (1)\quad y^{n+1} &= y^n+kAy^n,\\ (2)\quad y^{n+1} &= y^n+kAy^{n+1},\\ (3)\quad y^{n+1} &= y^n+\frac {k}{2}A\left (y^n+y^{n+1}\right ). \end {aligned} \]

Heat Equation Consider the following heat equation problem \[ u_t=u_{xx},\qquad t>0,\qquad 0<x<1, \] \[ u(0,t)=u(1,t)=0,\qquad t>0, \] \[ u(x,0)=u_0(x),\qquad 0\le x\le 1. \] Using the forward or backward difference for the time derivative and central difference for the space derivative, derive the explicit and implicit schemes. What’s the stability condition for each case? Perform a Fourier analysis on both the explicit and implicit methods and use the von Neumann condition to verify these stability conditions.

Find the radius of convergence of the Taylor series with the center at the origin (\(x_0=0\)) for the following function

1. \[ f(x)=\frac {\sqrt {x+3}}{x^2+3} \]
2. \[ f(x)=\frac {1}{(\cos x+3)^2} \]

(\(x\) is a real variable; you do not need to find the series themselves).

Define what is meant by a linear multistep method and its region of absolute stability. Show that a convergent explicit linear multistep method cannot have an infinite region of absolute stability.

Find the radius of convergence of the Taylor series with the center at \(x_0=1\) for the following function

1. \[ f(x)=x^{1/3} \]
2. \[ f(z)=\frac {(\sin x+2)^2}{(\sin x-2)^2} \]

(\(x\) is a real variable; you do not need to find the series themselves).

Find the radius of convergence of the Taylor series \[ \text {for the function } f(x)=\frac {1}{\cos ^2 x+4} \quad \text {around the origin } x_0=0. \] {You do not need to find the Taylor series itself.}

Heat Equation Stability:

Consider the variable coefficient diffusion equation \[ v_t=(\beta v_x)_x,\qquad 0<x<1,\quad t>0 \] with Dirichlet boundary conditions \[ v(0,t)=0,\qquad v(1,t)=0 \] and initial data \(v(x,0)=f(x)\). Assume that \(\beta (x)\geq \beta _0>0\), and that \(\beta (x)\) is smooth. Let \(\beta _{j+1/2}=\beta (x_{j+1/2})\). A scheme for this problem is: \[ \frac {u_j^{n+1}-u_j^n}{k} = \frac {1}{h^2} \left \{ \beta _{j-1/2}u_{j-1}^{n+1} -(\beta _{j-1/2}+\beta _{j+1/2})u_j^{n+1} +\beta _{j+1/2}u_{j+1}^{n+1} \right \}. \] Analyze the 2-norm stability of this scheme for solving this initial boundary value problem.

NOTE: Because of the boundary conditions and the variable coefficients, von Neumann (Fourier) analysis does not apply.

7) Numerical Methods for ODEs: Consider the initial value problem for a system of ODEs \[ y'=Ay \] where \(A\) has eigenvalues with real parts ranging from \(-10^6\) to \(-1\). Which of the following schemes would be the best choice for solving this problem? Justify your answer in terms of stability, accuracy, and efficiency. Does your answer change depending on the size of the system and the condition number of \(A\)? Explain your answer. \[ y^{n+1}=y^n+kAy^n, \tag {1} \] \[ y^{n+1}=y^n+kAy^{n+1}, \tag {2} \] \[ y^{n+1}=y^n+\frac {k}{2}\left (Ay^n+Ay^{n+1}\right ). \tag {3} \]

Find the radius of convergence of the Taylor series for the function \[ \frac {1}{2+\cos x} \qquad (x\text { is a real variable}) \] around the origin (\(x_0=0\)). {You do not need to find the Taylor series itself.}

Linear Multistep Methods

Construct an example of:

a) a consistent but not stable (not zero-stable) linear multistep method

b) a stable (zero-stable) but not consistent linear multistep method

What kind of behavior do you expect from the numerical solution produced by the methods in a) and in b)?

Convergence of Linear Multistep Method

Consider the method \[ y_{n+2}-2y_{n+1}+y_n = \frac {h}{2}\left (f(t_{n+2},y_{n+2})-f(t_n,y_n)\right ) \] Apply the method to the scalar IVP \(y'=y\), \(y(0)=1\) and solve exactly the resulting difference equation, considering the starting values to be \(y_0=y_1=1\). Show theoretically that the numerical solution does not converge as \(h\to 0\) and \(n\to \infty \).

Stability of the Scheme

Using the von Neumann method investigate the stability of the implicit downwind scheme: \[ \frac {u_m^{n+1}-u_m^n}{\Delta t} - \frac {u_m^{n+1}-u_{m-1}^{n+1}}{h} = f_m^n, \] \[ u_m^0=g_m,\qquad m=0,\pm 1,\pm 2,\ldots ,\qquad n=0,1,\ldots ,[T/\Delta t]-1. \]

Stability of the Scheme

Using the von Neumann method investigate the stability of the implicit upwind scheme: \[ \frac {u_m^{n+1}-u_m^n}{\Delta t} - \frac {u_{m+1}^{n+1}-u_m^{n+1}}{h} = f_m^n, \] \[ u_m^0=g_m,\qquad m=0,\pm 1,\pm 2,\ldots ,\qquad n=0,1,\ldots ,[T/\Delta t]-1. \] Comment on the CFL condition for this scheme.

Explicit Methods Define what we mean by Linear Multistep Method and what we mean by its region of absolute stability. Show that a convergent explicit linear multistep method cannot have an unbounded region of absolute stability.

Crank Nicolson Method Describe the Crank-Nicolson Method for the solution of the equation \[ u_t=u_{xx} \tag {6} \] and use the von Neumann method to examine its stability.

Local Truncation Error What is the local truncation error of a numerical scheme? What does consistency mean? Give a specific and precise example.

Iterative Methods for Linear Systems

Consider the boundary value problem \[ -u''(x)+\beta u(x)=f(x),\qquad \text {for } 0\leq x\leq 1 \] where \(\beta >0\), and with \(u(0)=u(1)=0\), and the following discretization of it: \[ -U_{j-1}+(2+\beta h^2)U_j-U_{j+1}=F_j \] for \(j=1,2,\ldots ,N-1\) where \(Nh=1\), \(F_j\equiv h^2f(jh)\), and \(U_0=U_N=0\).

Analyze the convergence properties of the Jacobi iterative method for this problem. In particular, express the speed of convergence as a function of the discretization stepsize \(h\). How does the number of iterations required to reduce the initial error by a factor \(\delta \) depend on \(h\)? In practice, would you use this method to solve the given problem? If so, explain why this is a good idea? If not, how would you solve it in practice?

Iterative Methods for Linear Systems

Consider the boundary value problem \[ -u''(x)+\beta u(x)=f(x),\qquad \text {for } 0\leq x\leq 1 \] where \(\beta >0\), and with \(u(0)=u(1)=0\), and the following discretization of it: \[ -U_{j-1}+(2+\beta h^2)U_j-U_{j+1}=F_j \] for \(j=1,2,\ldots ,N-1\) where \(Nh=1\), \(F_j\equiv h^2f(jh)\), and \(U_0=U_N=0\).

Analyze the convergence properties of the Jacobi iterative method for this problem. In particular, express the speed of convergence as a function of the discretization stepsize \(h\). How does the number of iterations required to reduce the initial error by a factor \(\delta \) depend on \(h\)? In practice, would you use this method to solve the given problem? If so, explain why this is a good idea? If not, how would you solve it in practice?

Absolute Stability Consider solving the initial value problem \[ y'=f(x,y),\qquad y(a)=y_0 \tag {4} \] by the linear multistep method (LMM) \[ \sum _{j=0}^k \alpha _j y_{n+j} = h\sum _{j=0}^k \beta _j f_{n+j}. \tag {5} \]

where we use the standard notation \[ x_n=a+nh,\qquad y_n\approx y(x_n),\qquad f_n=f(x_n,y_n). \tag {6} \] Define what it means that the LMM is explicit. Define its region of absolute stability. Show that an explicit linear multistep method cannot have an infinite region of absolute stability.

Convergence Consider the one-dimensional heat equation initial boundary value problem \[ u_t=u_{xx},\quad x\in [0,1],\quad t\geq 0,\quad u(x,0)=f(x),\quad u(0,t)=u(1,t)=0 \tag {7} \] and the discretization \[ x_m=mh,\qquad t_n=nk,\qquad \text {and}\qquad U_m^n\approx u(x_m,t_n) \tag {8} \] where \[ U_m^{n+1}=U_m^n+r\left (U_{m-1}^n-2U_m^n+U_{m+1}^n\right ) \tag {9} \] and \(r=k/h^2\) is the grid constant. Define what it means that this method is convergent, say for what values of \(r\) it is convergent, and show that your statement is correct.

ODEs

Consider the initial value problem for a system of ODEs \[ y'=Ay \] where \(A\) has eigenvalues with real parts ranging from \(-10^6\) to \(-1\). Which of the following schemes would be the best choice for solving this problem? Justify your answer in terms of stability, accuracy, and efficiency. Does your answer change depending on the size and condition number of \(A\)? Explain your answer. \[ y^{n+1}=y^n+kAy^n, \tag {1} \] \[ y^{n+1}=y^n+kAy^{n+1}, \tag {2} \] \[ y^{n+1}=y^n+\frac {k}{2}(Ay^n+Ay^{n+1}). \tag {3} \]

ODE Consider the initial value problem for a system of ODEs \[ y'=Ay \] where \(A\) has eigenvalues ranging from \(-10^4\) to \(-1\). Which of the following schemes would be the best choice for solving this problem? Justify your answer in terms of stability, accuracy, and efficiency. \[ y^{n+1}=y^n+hAy^n, \tag {1} \] \[ y^{n+1}=y^n+hAy^{n+1}, \tag {2} \] \[ y^{n+1}=y^n+\frac {h}{2}(Ay^n+Ay^{n+1}). \tag {3} \]

Elliptic Equations Consider the standard five-point difference approximation (centered difference for both the gradient and divergence operators) for the variable coefficient Poisson’s equation \[ -\nabla \cdot (a\nabla u)=f \] with Dirichlet boundary conditions, in a two-dimensional rectangular region. We assume that \(a(\mathbf {x})\ge \delta >0\). The approximate solution \(\{u_{i,j}\}\) satisfies a linear system \(Au=b\).

(a) State and prove the maximum principle for the numerical solution \(u_{i,j}\).

(b) Use the maximum principle to show that the coefficient matrix \(A\) is nonsingular.

(c) Derive the matrix \(A\) in the one-dimensional case and show that it is symmetric and positive definite.

(d) With this linear system, do we get convergence with Jacobi iterations? What about Gauss-Seidel iterations? Give reasons for your answers.

Wave Equation For the PDE \(v_t+cv_x=0\), derive the Lax-Wendroff scheme \[ u_j^{n+1}=u_j^n -\frac {1}{2}\left (\frac {c\Delta t}{\Delta x}\right ) (u_{j+1}^n-u_{j-1}^n) +\frac {1}{2}\left (\frac {c\Delta t}{\Delta x}\right )^2 (u_{j-1}^n-2u_j^n+u_{j+1}^n). \] Analyze this scheme in terms of accuracy and stability. Does this scheme satisfy a maximum principle? If it does, explain why. If it does not, give an example to show this.

Convection-Diffusion Equation For the convection-diffusion equation \[ v_t+av_x=bv_{xx}, \] where \(a>0\) and \(b>0\), a natural Crank-Nicolson scheme is \[ \begin {aligned} \frac {u_j^{n+1}-u_j^n}{\Delta t} &+\frac {a}{2}\left ( \frac {u_{j+1}^n-u_{j-1}^n}{2\Delta x} +\frac {u_{j+1}^{n+1}-u_{j-1}^{n+1}}{2\Delta x} \right ) \\ &= \frac {b}{2}\left ( \frac {u_{j+1}^n-2u_j^n+u_{j-1}^n}{(\Delta x)^2} +\frac {u_{j+1}^{n+1}-2u_j^{n+1}+u_{j-1}^{n+1}}{(\Delta x)^2} \right ). \end {aligned} \] Let \(\nu =b\Delta t/(\Delta x)^2\) and \(\mu =a\Delta t/(2\Delta x)\). Show that the scheme is stable in the maximum norm if \(\mu \le \nu \le 1\). Also use Fourier analysis to derive the stability condition in the \(\ell _2\) norm and show that the scheme is always stable in the \(\ell _2\) norm. Suppose \(u\) stands for some volume fraction of a certain species, which norm is more appropriate? Explain why this scheme is so restrictive that it is not a good choice for many applications.

Ordinary Differential Equations Consider the initial value problem for a system of two ODEs \[ y'=Ay \] where \(A\) has eigenvalues \(-10^4\) and \(1\). Which of the following schemes would be the best choice for solving this problem? Justify your answer in terms of stability, accuracy, and efficiency. \[ y^{n+1}=y^n+hAy^n, \tag {1} \] \[ y^{n+1}=y^n+hAy^{n+1}, \tag {2} \] \[ y^{n+1}=y^n+\frac {h}{2}(Ay^n+Ay^{n+1}). \tag {3} \]

Wave Equation Consider the Lax-Friedrichs scheme below for the PDE \(v_t+cv_x=0\), where you may assume that \(c>0\): \[ u_j^{n+1}=\frac {1}{2}(u_{j-1}^n+u_{j+1}^n) -\frac {ck}{2h}(u_{j+1}^n-u_{j-1}^n), \] where \(h\) is the spacestep and \(k\) is the timestep. Note that the scheme can also be written \[ u_j^{n+1}=u_j^n-\frac {ck}{2h}(u_{j+1}^n-u_{j-1}^n) +\frac {1}{2}(u_{j-1}^n-2u_j^n+u_{j+1}^n). \]

a) Analyze the stability of this scheme in terms of \(\alpha =\frac {ck}{h}\).

b) Recall that a modified equation for a scheme is the PDE satisfied, to a certain order, by a smooth solution \(u(x,t)\) of the difference scheme. To zeroth order in \(k\) and \(h\), the modified equation for the Lax-Friedrichs scheme is \(u_t+cu_x=0\). Derive the modified equation through terms of first order in \(k\) and \(h\). You may assume that \(\alpha \) is constant. Does this modified equation predict the same stability condition as your analysis in part a)? Explain your answer. What does the modified equation suggest about the numerical solution that will be obtained for the step-function initial data \(u(x,0)=1\) for \(x\le 0\) and \(u(x,0)=0\) for \(x>0\)?

Heat Equation Consider the variable coefficient diffusion equation \[ v_t=(\beta v_x)_x,\qquad 0<x<1,\quad t>0 \] with Dirichlet boundary conditions \[ v(0,t)=0.0,\qquad v(1,t)=0.0 \] and initial data \(v(x,0)=f(x)\). Assume that \(\beta (x)\ge \beta _0>0\), and that \(\beta (x)\) is smooth. For any fixed value of \(\theta \in [0,1]\), consider the \(\theta \)-scheme for this problem: \[ \begin {aligned} \frac {u_j^{n+1}-u_j^n}{k} =\frac {1}{h^2}\{& (1-\theta )(\beta _{j-1/2}u_{j-1}^{n+1} -(\beta _{j-1/2}+\beta _{j+1/2})u_j^{n+1} +\beta _{j+1/2}u_{j+1}^{n+1})\\ &+\theta (\beta _{j-1/2}u_{j-1}^n -(\beta _{j-1/2}+\beta _{j+1/2})u_j^n +\beta _{j+1/2}u_{j+1}^n)\}. \end {aligned} \] with appropriate treatment of the boundary conditions in the relevant equations. Here, \(h\) is the spacestep, \((J+1)h=1\), \(k\) is the timestep, and \(\beta _{j+1/2}=\beta (x_{j+1/2})\), where \(x_{j+1/2}=(j+1/2)h\). Let \(M\) be the symmetric \(J\times J\) matrix defined by \[ \setlength {\arraycolsep }{3pt} M= \begin {pmatrix} -(\beta _{1/2}+\beta _{3/2}) & \beta _{3/2} & \cdots & \cdots & \cdots & \cdots \\ \beta _{3/2} & -(\beta _{3/2}+\beta _{5/2}) & \beta _{5/2} & \cdots & \cdots & \cdots \\ \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \\ \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \\ \cdots & \cdots & \cdots & \beta _{J-3/2} & -(\beta _{J-3/2}+\beta _{J-1/2}) & \beta _{J-1/2}\\ \cdots & \cdots & \cdots & \cdots & \beta _{J-1/2} & -(\beta _{J-1/2}+\beta _{J+1/2}) \end {pmatrix}. \] Analyze the stability of this scheme in terms of \(\theta \), \(\frac {k}{h^2}\), and the eigenvalues of \(M\). Are there values of \(\theta \) for which the scheme is stable for all choices of \(\frac {k}{h^2}>0\)? (Note that answering this question should not require lots of writing.)

Eigenvalue Problems

a) Let \(A\) be a real \(n\times n\) matrix with simple eigenvalue \(\lambda \), right eigenvector \(x\) (i.e., \(Ax=\lambda x\)), and left eigenvector \(y\) (i.e., \(y^*A=\lambda y^*\)). Assume that \(\lVert x\rVert _2=\lVert y\rVert _2=1\). Show that the condition number for \(\lambda \) to changes in \(A\) is given by \(s(\lambda )=1/\lvert y^*x\rvert \). Give an example of a matrix \(A\) which has a badly conditioned eigenvalue, and another matrix \(A\) which has only well-conditioned eigenvalues.

b) Consider the problem of computing an eigenvalue/eigenvector pair of a real symmetric \(m\times m\) matrix \(A\) with eigenvalues of distinct magnitudes. One algorithm for trying to do this is Rayleigh Quotient Iteration:

Guess \(v^{(0)}\) with \(\lVert v^{(0)}\rVert _2=1\).

Set \(\lambda ^{(0)}=(v^{(0)})^TAv^{(0)}\).

For \(k=1,2,3,\ldots \): \[ \text {Solve }(A-\lambda ^{(k-1)}I)w=v^{(k-1)}\text { for }w, \] \[ v^{(k)}=\frac {w}{\lVert w\rVert _2}, \] \[ \lambda ^{(k)}=(v^{(k)})^TAv^{(k)}. \] Until convergence.

Suppose that the Rayleigh Quotient Iteration iterates \(\lambda ^{(k)}\) converge to a simple eigenvalue \(\lambda \) of \(A\). Analyze the speed of convergence of \(\lambda ^{(k)}\) to \(\lambda \).

Eigenvalue Problems

a) Let \(A\) be a real \(n\times n\) matrix with simple eigenvalue \(\lambda \), right eigenvector \(x\) (i.e., \(Ax=\lambda x\)), and left eigenvector \(y\) (i.e., \(y^*A=\lambda y^*\)). Assume that \(\lVert x\rVert _2=\lVert y\rVert _2=1\). Show that the condition number for \(\lambda \) to changes in \(A\) is given by \(s(\lambda )=1/\lvert y^*x\rvert \). Give an example of a matrix \(A\) which has a badly conditioned eigenvalue, and another matrix \(A\) which has only well-conditioned eigenvalues.

b) Consider the problem of computing an eigenvalue/eigenvector pair of a real symmetric \(m\times m\) matrix \(A\) with eigenvalues of distinct magnitudes. One algorithm for trying to do this is Rayleigh Quotient Iteration:

Guess \(v^{(0)}\) with \(\lVert v^{(0)}\rVert _2=1\).

Set \(\lambda ^{(0)}=(v^{(0)})^TAv^{(0)}\).

For \(k=1,2,3,\ldots \): \[ \text {Solve }(A-\lambda ^{(k-1)}I)w=v^{(k-1)}\text { for }w, \] \[ v^{(k)}=\frac {w}{\lVert w\rVert _2}, \] \[ \lambda ^{(k)}=(v^{(k)})^TAv^{(k)}. \] End.

Suppose that the Rayleigh Quotient Iteration iterates \(\lambda ^{(k)}\) converge to a simple eigenvalue \(\lambda \) of \(A\). Analyze the speed of convergence of \(\lambda ^{(k)}\) to \(\lambda \).

Iterative Methods for Linear Systems Consider the system \(Ax=b\), in which \(A\) is a symmetric positive definite real \(m\times m\) matrix, and denote by \(x_*\) its true solution. The conjugate gradient (CG) algorithm for finding the solution to the system is:

Set \(x_0=0\), \(r_0=b\), and \(p_0=r_0\).

For \(n=1,2,3,\ldots \), \[ \alpha _n=\frac {r_{n-1}^T r_{n-1}}{p_{n-1}^T A p_{n-1}}, \] \[ x_n=x_{n-1}+\alpha _n p_{n-1}, \] \[ r_n=r_{n-1}-\alpha _n A p_{n-1}, \] \[ \beta _n=\frac {r_n^T r_n}{r_{n-1}^T r_{n-1}}, \] \[ p_n=r_n+\beta _n p_{n-1}. \] End.

The vectors generated by this algorithm have the property that, provided \(r_{n-1}\ne 0\), \[ \begin {aligned} \mathcal {K}_n &= \operatorname {span}\{b,Ab,A^2b,\ldots ,A^{n-1}b\} = \operatorname {span}\{x_1,x_2,\ldots ,x_n\}\\ &= \operatorname {span}\{r_0,r_1,\ldots ,r_{n-1}\} = \operatorname {span}\{p_0,p_1,\ldots ,p_{n-1}\} \end {aligned} \] and \(r_n^T r_j=0\) and \(p_n^T A p_j=0\) for all \(j<n\).

Let \(e_n=\lVert x_n-x_*\rVert _A\), where the \(A\)-norm of a vector \(x\) is defined by \(\lVert x\rVert _A=(x^T Ax)^{1/2}\).

a) Prove the following: Suppose the CG algorithm is applied to the symmetric positive definite system \(Ax=b\) and that it has not yet converged (\(r_{n-1}\ne 0\)). Then \(e_n<\lVert y-x_*\rVert _A\) for all \(y\in \mathcal {K}_n\), \(y\ne x_n\); \(\lVert e_n\rVert _A\le \lVert e_{n-1}\rVert _A\); and \(e_N\to 0\) for some \(N\le m\).

b) In practice this result is not very useful. Why is this, and what determines the speed of convergence of \(x_n\) to \(x_*\)? What is ‘preconditioning’ and how is it related to the speed of convergence of the CG algorithm?

Eigenvalue Problems

a) (8%) QR iteration can be viewed as an extension of the power method. We assume \(n\times n\) real symmetric, nonsingular matrix \(A\), with distinct eigenvalues. Here is the algorithm of orthogonal iteration for two linearly independent vectors \(z_1\) and \(z_2\), listed in an \(n\times 2\) matrix \(Z_0=[z_1,z_2]\):

\(i=0\)

repeat \[ \begin {aligned} Y_{i+1} &= AZ_i,\\ \text {QR Factor }Y_{i+1} &= Z_{i+1}R_{i+1},\\ i &= i+1 \end {aligned} \] until convergence

Notice that \(Y_i,Z_i\) are \(n\times 2\) matrices, and \(R_i\) is a \(2\times 2\) upper triangular matrix. Our assumptions allow us to write \(A=S\Lambda S^{-1}\), where \(\Lambda \) is a diagonal matrix containing all eigenvalues \(\lvert \lambda _1\rvert >\lvert \lambda _2\rvert >\cdots >\lvert \lambda _n\rvert \), \(S=[s_1,s_2,\ldots ,s_n]\) contains all the eigenvectors.

Show that \[ \operatorname {span}\{Z_{i+1}\} = \operatorname {span}\{Y_{i+1}\} = \operatorname {span}\{AZ_i\}, \] so \[ \operatorname {span}\{Z_i\} = \operatorname {span}\{A^iZ_0\} = \operatorname {span}\{S\Lambda ^iS^{-1}Z_0\}. \] Analyze it by writing \[ S\Lambda ^iS^{-1}Z_0=\lambda _2^i S D_i(S^{-1}Z_0). \] Derive the diagonal \(D_i\) and the \(n\times 2\) matrix \(Z'_0=S^{-1}Z_0\), find the limiting span of \(Z_i\) as \(i\to \infty \). How do we obtain the approximating eigenvalues \(\lambda _1,\lambda _2\) and approximating eigenvectors \(s_1\) and \(s_2\), once iterations are terminated?

b) (6%) Assume that we can continue this extension to \(n\) linearly independent vectors, contained in an \(n\times n\) nonsingular \(Z_0\); in particular when \(Z_0=I\), we will have the method of QR iteration. Describe the QR algorithm (without shift) and argue for the rate of convergence based on the analysis in part a).

c) (6%) To speed up convergence, we should use shifting. Describe the QR iteration with a shift, and make a suggestion as what shift should be used. How fast is the convergence?

Numerical ODEs

a) (8%) Consider the general system of ODEs for \(y(x)\in \mathbb {R}^n\), \(x\in \mathbb {R}\), and \[ f=[f_1(x,y),\ldots ,f_n(x,y)]^T: \qquad y'=f(x,y), \] where we assume \(f\) is continuously differentiable in both \(x\) and \(y\). To study local stability, we compare two solutions \(y_1(x)\) and \(y_2(x)\) where \(y_1(x_0)\approx y_2(x_0)\). Derive the linearized system for \(e(x)=y_1(x)-y_2(x)\) near \(x=x_0\) to show that \[ e'\approx Ae \] and find this square matrix \(A\). What is the most important property of \(A\) we need to find that will affect the choice of numerical schemes and step sizes?

b) (12%) Now consider the initial value problem \[ y''+2\mu y'+(1+\mu ^2)y=0, \] \[ y(0)=0,\qquad y'(0)=1, \] where \(\mu >0\) is a given constant. If we use a simple explicit Euler’s method, determine the maximum step size \(h\) that is allowed for the solution to be stable.

Wave Equations Consider the standard wave equation in one space dimension: \[ u_{tt}=c^2u_{xx},\qquad 0<x<1,\quad t>0, \] with initial conditions \[ u(x,0)=f(x),\qquad u_t(x,0)=g(x),\qquad 0<x<1 \] for given smooth functions \(f\) and \(g\), and boundary condition \[ u(0,t)=u(1,t)=0,\qquad t>0. \] The leap-frog scheme takes advantage of the fact that the equation can be reformulated as a system of first order equations \[ u_t+cv_x=0, \qquad v_t+cu_x=0. \] Denote \(x_j=j\Delta x\), \(U_j^n\approx u(x_j,n\Delta t)\), and \(V_{j+1/2}^{n+1/2}\approx v(x_{j+1/2},(n+\frac {1}{2})\Delta t)\); use central difference approximations for both derivatives to derive the leap-frog scheme for the system of equations. How do you implement the initial and boundary conditions for both \(u\) and \(v\)? Derive the truncation error for this approximation and determine the order of truncation error. For stability, consider solutions in the form \[ \begin {bmatrix} U^n\\ V^{n+1/2} \end {bmatrix} = \lambda ^n e^{ikx} \begin {bmatrix} \widetilde {U}\\ \widetilde {V} \end {bmatrix}, \] for some constants \(\widetilde {U}\) and \(\widetilde {V}\). Substitute this solution into the finite difference equations, find out what condition \(\lambda \) has to satisfy in order for \(\widetilde {U}\) and \(\widetilde {V}\) to have nontrivial solutions. For those \(\lambda \) values, compute \(\lvert \lambda \rvert \) to see if the stability condition \(\lvert \lambda \rvert \le 1\) is satisfied.

Convection-Diffusion Equation For the convection-diffusion equation \[ u_t+au_x=bu_{xx}, \] where \(a>0\) and \(b>0\), a natural Crank-Nicolson scheme is \[ \begin {aligned} \frac {U_j^{n+1}-U_j^n}{\Delta t} &+\frac {a}{2}\left ( \frac {U_{j+1}^n-U_{j-1}^n}{2\Delta x} + \frac {U_{j+1}^{n+1}-U_{j-1}^{n+1}}{2\Delta x} \right )\\ &= \frac {b}{2}\left ( \frac {U_{j+1}^n-2U_j^n+U_{j-1}^n}{(\Delta x)^2} + \frac {U_{j+1}^{n+1}-2U_j^{n+1}+U_{j-1}^{n+1}}{(\Delta x)^2} \right ). \end {aligned} \] Denote \(\nu =b\Delta t/(\Delta x)^2\) and \(\mu =a\Delta t/(2\Delta x)\); show that the scheme is stable in the max norm if \(\mu \le \nu \le 1\). Also use Fourier analysis to derive the stability condition in the \(\ell _2\) norm and show that the scheme is always stable in the \(\ell _2\) norm. Suppose \(u\) stands for some chemical concentration, which norm is more appropriate? Explain why this scheme is so restrictive that it is not often used.

ODEs Consider the initial value problem of ODEs: \[ y'=f(x,y),\qquad y(a)=y_0. \tag {10} \] Describe the Trapezoidal Rule applied to this problem and analyze its local accuracy and its region of absolute stability.

PDEs Consider the one-dimensional Heat Equation \[ u_t=u_{xx},\qquad t\ge 0,\qquad x\in [0,1],\qquad u(x,0)=f(x),\qquad u(0,t)=u(1,t)=0. \tag {11} \] Describe how to discretize this problem by applying the Method of Lines and the Backward Euler Method. Analyze the stability and local accuracy of your method.

Numerical ODEs Define the terms linear multistep method (LMM), order, consistency, zero-stability, local truncation error, and region of absolute stability of an LMM. What is the highest order implicit 1 step method and discuss its region of absolute stability?

Numerical PDEs Consider the Initial-Boundary Value Problem \[ u_t=u_{xx},\qquad x\in [0,1],\qquad t\ge 0, \qquad u\text { given for }x=0,\ x=1,\text { and }t=0. \tag {4} \] Describe the Crank-Nicolson Method, and analyze its local truncation error and stability.

Eigenvalue Problems

(a) For the general eigenvalue problem of a nonsingular matrix \(A\), the simplest method is the power method, which converges to the largest eigenvalue in magnitude, assuming it is simple. The algorithm is the following:

\(i=0\)

repeat \[ \begin {aligned} y_{i+1} &= Ax_i,\\ x_{i+1} &= y_{i+1}/\lVert y_{i+1}\rVert _2,\\ \tilde {\lambda }_{i+1} &= x_{i+1}^T A x_{i+1},\\ i &= i+1 \end {aligned} \] until convergence

The inverse iteration algorithm takes advantage of the fact the if \(\sigma \) (the shift) is closest to one of the simple eigenvalues \(\lambda \), then \(1/(\lambda -\sigma )\) will be the largest eigenvalue in magnitude of the matrix \((A-\sigma I)^{-1}\). Use this fact to construct the inverse iteration algorithm. Note that we will not need to construct the explicit inverse of \(A-\sigma I\). How fast does it converge to the desired eigenvalue? When do you expect the convergence to be slow? How do we find the corresponding eigenvector?

(b) Now consider the symmetric eigenvalue problem where \(A\) is symmetric. What is the Rayleigh quotient? The Rayleigh quotient method is the inverse iteration algorithm described above applied to a symmetric matrix, when the shift is chosen to be the Rayleigh quotient \(\rho _i\) of \(A\) and the updated vector \(x_i\). Explain the rationale of choosing \(\rho _i\) to be the shift and determine the speed of convergence of \(\rho _i\) to \(\lambda \).

Numerical ODEs

(a) Consider the general system of ODEs for \(y(x)\in \mathbb {R}^n\), \(x\in \mathbb {R}\) and \(f=[f_1,\ldots ,f_n]^T\): \[ y'=f(x,y), \] where we assume \(f\) to be Lipschitz continuous. To study the numerical stability, we linearize \(f\) near \(x=x_0,y=y_0\) to obtain: \[ y'\approx f(x_0,y_0)+f_x(x_0,y_0)(x-x_0)+\nabla _y f(x_0,y_0)(y-y_0). \] If we denote \(A=\nabla _y f(x_0,y_0)\), then we claim that the local stability of the above system can be analyzed by studying the following constant coefficient system \[ y'=Ay. \] Justify this claim by analyzing the simple scalar equation \[ y'=a+bx+cy, \] where \(a,b,c\) are constants, in the explicit Euler’s method, to show that the stability is determined by how we solve \(y'=cy\). What is the most important property of \(A\) we need to find that will affect the choice of numerical schemes and step sizes?

(b) Now consider the initial value problem \[ y''+2\mu y'+(1+\mu ^2)y=0, \] \[ y(0)=0,\qquad y'(0)=1, \] where \(\mu >0\) is a given constant. If we use a simple explicit Euler’s method, determine the maximum step size \(h\) that is allowed for the solution to be stable.

Wave Equations Consider the first-order wave equation in one space dimension: \[ u_t+au_x=0,\qquad -\infty <x<\infty ,\qquad t>0, \] where \(a\) is a constant, with the initial condition \[ u(x,0)=u_0(x),\qquad -\infty <x<\infty \] for a given function \(u_0\). We use the usual discretization \[ U_j^n\approx u(x_j,t^n), \] where \(x_j=j\Delta x\) and \(t^n=n\Delta t\). For a finite difference scheme to solve this problem, what is the CFL condition? Give an intuitive explanation for the condition. Derive the first-order upwind scheme for this equation and its truncation error. For numerical stability, derive the conditions for stability in the maximum norm (maximum principle) and the \(\ell _2\) norm (Fourier analysis).

Heat Equation Derive the Crank-Nicolson scheme for the following problem \[ u_t=u_{xx},\qquad t>0,\qquad 0<x<1, \] with the boundary conditions \[ u(0,t)=u(1,t)=0,\qquad t\ge 0, \] and the initial data \[ u(x,0)=u_0(x),\qquad 0\le x\le 1, \] for some given function \(u_0\). Use the usual discretization with uniform mesh points \(x_j=jh\) where \(h=1/N\) and \(j=0,1,2,\ldots ,N\). Show that the method is second order by deriving the truncation error. What other conditions do we need to guarantee the convergence of the numerical solution? Derive the stability condition in \(\ell _2\) norm by using Fourier analysis. For any value of \(\nu =\Delta t/(\Delta x)^2\), do we always have the maximum principle satisfied? If we have some initial data such that \(0\le u_0(x)\le 1\), derive the condition that should be imposed on \(\nu \) so that \(0\le u_j^n\le 1\) also holds for all \(j\) and \(n>0\).

Iterative Methods for Linear Systems Consider the fixed-point iteration \[ u^{(k+1)}=Tu^{(k)}+c \] for finding a solution of the problem \[ u=Tu+c \] where \(T\) is an \(m\times m\) real matrix and \(c\) is a real \(m\)-vector.

a) Show that the fixed point iteration will converge for an arbitrary initial guess \(u^{(0)}\) if and only if the spectral radius of \(T\), \(\rho (T)\), is less than 1.

b) Consider the boundary value problem \[ -u''(x)=f(x),\qquad \text {for }0\le x\le 1 \] with \(u(0)=u(1)=0\), and the following discretization of it: \[ -U_{j-1}+2U_j-U_{j+1}=F_j \] for \(j=1,2,\ldots ,N-1\) where \(Nh=1\) and \(F_j=h^2f(jh)\).

Analyze the convergence properties of the Jacobi iterative method for this problem. In particular, express the speed of convergence as a function of the discretization stepsize \(h\). How does the number of iterations required to reduce the initial error by a factor of \(\delta \) depend on \(h\)? In practice, would you use this method to solve the given problem? If so, explain why this is a good idea? If not, how would you solve it in practice?

Eigenvalue Problems

a) Let \(A\) be a real \(n\times n\) matrix with simple eigenvalue \(\lambda \), right eigenvector \(x\) (i.e., \(Ax=\lambda x\)), and left eigenvector \(y\) (i.e., \(y^*A=\lambda y^*\)). Assume that \(\lVert x\rVert _2=\lVert y\rVert _2=1\). Show that the condition number for \(\lambda \) to changes in \(A\) is given by \(s(\lambda )=1/\lvert y^*x\rvert \). Give an example of a matrix \(A\) which has a badly conditioned eigenvalue, and another matrix \(A\) which has only well-conditioned eigenvalues.

b) Consider the problem of computing an eigenvalue/eigenvector pair of a real symmetric \(m\times m\) matrix \(A\) with eigenvalues of distinct magnitudes. One algorithm for trying to do this is Rayleigh Quotient Iteration:

Guess \(v^{(0)}\) with \(\lVert v^{(0)}\rVert _2=1\).

Set \(\lambda ^{(0)}=(v^{(0)})^TAv^{(0)}\).

For \(k=1,2,3,\ldots \) \[ \begin {aligned} &\text {Solve }(A-\lambda ^{(k-1)}I)w=v^{(k-1)}\text { for }w,\\ &v^{(k)}=w/\lVert w\rVert _2,\\ &\lambda ^{(k)}=(v^{(k)})^TAv^{(k)}. \end {aligned} \] End

Suppose that the Rayleigh Quotient Iteration iterates \(\lambda ^{(k)}\) converge to a simple eigenvalue \(\lambda \) of \(A\). Analyze the speed of convergence of \(\lambda ^{(k)}\) to \(\lambda \).

Krylov Subspace Methods Consider the system \(Ax=b\), in which \(A\) is a symmetric positive definite real \(m\times m\) matrix, and denote by \(x_*\) its true solution. The conjugate gradient (CG) algorithm for finding the solution to the system is:

Set \(x_0=0\), \(r_0=b\), and \(p_0=r_0\).

For \(n=1,2,3,\ldots \) \[ \begin {aligned} \alpha _n &= r_{n-1}^T r_{n-1}/p_{n-1}^T A p_{n-1},\\ x_n &= x_{n-1}+\alpha _n p_{n-1},\\ r_n &= r_{n-1}-\alpha _n A p_{n-1},\\ \beta _n &= r_n^T r_n/r_{n-1}^T r_{n-1},\\ p_n &= r_{n-1}+\beta _n p_{n-1}. \end {aligned} \] End

The vectors generated by this algorithm have the property that, provided \(r_{n-1}\ne 0\), \[ \begin {aligned} \mathcal {K}_n &\equiv \operatorname {span}\{b,Ab,A^2b,\ldots ,A^{n-1}b\} =\operatorname {span}\{x_1,x_2,\ldots ,x_n\}\\ &=\operatorname {span}\{r_0,r_1,\ldots ,r_{n-1}\} =\operatorname {span}\{p_0,p_1,\ldots ,p_{n-1}\} \end {aligned} \] and \(r_n^T r_j=0\) and \(p_n^T A p_j=0\) for all \(j<n\).

Let \(e_n=\lVert x_n-x_*\rVert _A\), where the \(A\)-norm of a vector \(x\) is defined by \(\lVert x\rVert _A=(x^TAx)^{1/2}\).

a) Prove the following:

Suppose the CG algorithm is applied to the symmetric positive definite system \(Ax=b\) and that it has not yet converged (\(r_{n-1}\ne 0\)). Then \(e_n<\lVert y-x_*\rVert _A\) for all \(y\in \mathcal {K}_n\), \(y\ne x_n\); \(\lVert e_n\rVert _A\le \lVert e_{n-1}\rVert _A\); and \(e_N=0\) for some \(N\le m\).

b) In practice this result is not very useful. Why is this, and what determines the speed of convergence of \(x_n\) to \(x_*\)? What is ‘preconditioning’ and how is it related to the speed of convergence of the CG algorithm?

Numerical Methods for ODEs Consider the Linear Multistep Method \[ y_{n+2}-\frac {4}{3}y_{n+1}+\frac {1}{3}y_n=\frac {2}{3}hf_{n+2} \] for solving an initial value problem \(y'=f(y,x)\), \(y(0)=\eta \). You may assume that \(f\) is Lipschitz continuous with respect to \(y\) for all \(x\).

a) Analyze the consistency, stability, accuracy, and convergence properties of this method.

b) Would it be more reasonable to use this method or Euler’s method for the initial value problem: \[ y'=-10^8[y-\cos (x)]-\sin (x),\qquad y(0)=1? \] Justify your answer.

Iterative Methods for Linear Systems Consider the linear system \(Ax=b\) and a splitting of \(A:A=M-K\), with \(M\) nonsingular, and the iteration \[ x_{m+1}=Rx_m+c, \] where \(R=M^{-1}K\) and \(c=M^{-1}b\). Show that the method converges for any \(x_0\) if \(\lVert R\rVert <1\) where \(\lVert \cdot \rVert \) is any operator norm.

Given an \((N-1)\times (N-1)\) \[ \setlength {\arraycolsep }{5pt} A= \begin {bmatrix} 2 & -1 & 0 & \cdots & 0\\ -1 & 2 & -1 & \ddots & 0\\ 0 & -1 & \ddots & \ddots & 0\\ \vdots & \ddots & -1 & 2 & -1\\ 0 & \cdots & 0 & -1 & 2 \end {bmatrix}, \qquad b=\frac {1}{N^2} \begin {pmatrix} f_1\\ f_2\\ \cdot \\ f_{N-2}\\ f_{N-1} \end {pmatrix}. \] Analyze the convergence properties of the Jacobi and Gauss-Seidel methods for this problem. For each method, express the convergence speed as a function of \(N\). How does the number of iterations required to reduce the initial error by a factor of \(\delta \) depend on \(N\) for each method?

Heat Equation Consider the following heat equation problem \[ u_t=u_{xx},\qquad t>0,\qquad 0<x<1, \] \[ u(0,t)=u(1,t)=0,\qquad t>0, \] \[ u(x,0)=u_0(x),\qquad 0\le x\le 1. \] Using the forward or backward difference for the time derivative and central difference for the space derivative, derive the explicit and implicit schemes. What’s the stability condition for each case? Perform a Fourier analysis on both the explicit and implicit methods and use the von Neumann condition to verify these stability conditions.

Numerical ODEs Define the terms linear multistep method (LMM), order, consistency, zero-stability, local truncation error, and region of absolute stability of an LMM. Derive the highest order implicit 2 step method and discuss its region of absolute stability.

Numerical PDEs Consider the PDE \[ u_t+u_x=0 \tag {5} \] and the usual discretization \[ U_m^n\approx u(x_m,t_n) \qquad \text {where}\qquad x_m=mh\quad \text {and}\quad t_n=nk. \tag {6} \] The Leap Frog Method is given by \[ U_m^{n+1}=U_m^{n-1}-\frac {k}{h}\left (U_{m+1}^n-U_{m-1}^n\right ). \tag {7} \] Compute the local truncation error and discuss the stability of this method.

Iterative Methods for Linear Systems a) Consider the fixed-point iteration \[ u^{(k+1)}=Tu^{(k)}+c \] for finding a solution of the problem \[ u=Tu+c \] where \(T\) is an \(m\times m\) real matrix and \(c\) is a real \(m\)-vector. Show that the fixed point iteration will converge for an arbitrary initial guess \(u^{(0)}\) if and only if the spectral radius of \(T\), \(\rho (T)\), is less than 1.

b) Consider the boundary value problem \[ -u''(x)=f(x),\qquad \text {for }0\le x\le 1 \] with \(u(0)=u(1)=0\), and the following discretization of it: \[ -U_{j-1}+2U_j-U_{j+1}=F_j \] for \(j=1,2,\ldots ,N-1\) where \(Nh=1\) and \(F_j=h^2f(jh)\).

Analyze the convergence properties of the Jacobi iterative method for this problem. In particular, express the speed of convergence as a function of the discretization stepsize \(h\). How does the number of iterations required to reduce the initial error by a factor of \(\delta \) depend on \(h\)? In practice, would you use this method to solve the given problem? If so, explain why this is a good idea? If not, how would you solve it in practice?

Eigenvalue Problems

a) Consider the problem of computing an eigenvalue/eigenvector pair of a real symmetric \(m\times m\) matrix \(A\) with eigenvalues of distinct magnitudes. Two algorithms for trying to do this are Power Iteration and Rayleigh Quotient Iteration:

i) Analyze the convergence properties of the Power Iteration algorithm for a generic guess \(v^0\).

ii) Suppose that the Rayleigh Quotient Iteration iterates \(\lambda ^{(k)}\) converge to a simple eigenvalue \(\lambda \) of \(A\). Analyze the speed of convergence of \(\lambda ^{(k)}\) to \(\lambda \).

b) Many algorithms for finding eigenvalues of a complex \(m\times m\) matrix \(A\) are based on the existence of a Schur decomposition of \(A\), namely, \(A=QTQ^*\) where \(Q\) is unitary and \(T\) is upper triangular. Prove that every matrix \(A\) has a Schur decomposition and explain how knowing the Schur decomposition of \(A\) would help one determine the eigenvalues of \(A\).

Krylov Subspace Methods Consider the system \(Ax=b\), in which \(A\) is a symmetric positive definite real \(m\times m\) matrix, and denote by \(x_*\) its true solution. The conjugate gradient (CG) algorithm for finding the solution to the system is:

Set \(x_0=0\), \(r_0=b\), and \(p_0=r_0\).

For \(n=1,2,3,\ldots \) \[ \begin {aligned} \alpha _n &= r_{n-1}^T r_{n-1}/p_{n-1}^T A p_{n-1},\\ x_n &= x_{n-1}+\alpha _n p_{n-1},\\ r_n &= r_{n-1}-\alpha _n A p_{n-1},\\ \beta _n &= r_n^T r_n/r_{n-1}^T r_{n-1},\\ p_n &= r_{n-1}+\beta _n p_{n-1}. \end {aligned} \] End

The vectors generated by this algorithm have the property that, provided \(r_{n-1}\ne 0\), \[ \begin {aligned} \mathcal {K}_n &\equiv \operatorname {span}\{b,Ab,A^2b,\ldots ,A^{n-1}b\} =\operatorname {span}\{x_1,x_2,\ldots ,x_n\}\\ &=\operatorname {span}\{r_0,r_1,\ldots ,r_{n-1}\} =\operatorname {span}\{p_0,p_1,\ldots ,p_{n-1}\} \end {aligned} \] and \(r_n^T r_j=0\) and \(p_n^T A p_j=0\) for all \(j<n\).

Let \(e_n=\lVert x_n-x_*\rVert _A\), where the \(A\)-norm of a vector \(x\) is defined by \(\lVert x\rVert _A=(x^TAx)^{1/2}\).

a) Prove the following:

Theorem: Suppose the CG algorithm is applied to the symmetric positive definite system \(Ax=b\) and that it has not yet converged (\(r_{n-1}\ne 0\)). Then \(e_n<\lVert y-x_*\rVert _A\) for all \(y\in \mathcal {K}_n\), \(y\ne x_n\); \(\lVert e_n\rVert _A\le \lVert e_{n-1}\rVert _A\); and \(e_N=0\) for some \(N\le m\).

b) In practice this result is not very useful. Why is this, and what determines the speed of convergence of \(x_n\) to \(x_*\)? What is ‘preconditioning’ and how is it related to the speed of convergence of the CG algorithm?

Numerical Methods for ODEs Consider the Linear Multistep Method \[ y_{n+1}-y_n=\frac {1}{2}h(f_{n+1}+f_n) \] for solving an initial value problem \(y'=f(y,x)\), \(y(0)=\eta \). You may assume that \(f\) is Lipschitz continuous with respect to \(y\) for all \(x\).

a) Analyze the consistency, stability, accuracy, and convergence properties of this method.

b) Would it be more reasonable to use this method or Euler’s method for the initial value problem: \[ y'=-10^8[y-\cos (x)]-\sin (x),\qquad y(0)=1? \] Justify your answer.

Iterative Methods for Linear Systems Consider the boundary value problem \[ -u''(x)=f(x),\qquad \text {for }0\le x\le 1 \] with \(u(0)=u(1)=0\), and the following discretization of it: \[ -U_{j-1}+2U_j-U_{j+1}=F_j \] for \(j=1,2,\ldots ,N-1\) where \(Nh=1\) and \(F_j=h^2f(jh)\).

Analyze the convergence properties of both the Jacobi and Gauss Seidel iterative methods for this problem. In particular, for each, express the speed of convergence as a function of the discretization stepsize \(h\). How does the number of iterations required to reduce the initial error by a factor \(\delta \) depend on \(h\) for each method? In practice, would you use one of these methods to solve the given problem? If so, explain why this is a good idea? If not, how would you solve it in practice?

Newton’s Method Suppose \(f(x)\) is a real-valued function defined for all values of the real variable \(x\), and that \(f'\) is Lipschitz continuous. State and prove a local convergence theorem for Newton’s method for finding a simple zero of \(f\). How does this result change if the root is a double root? Is Newton’s method globally convergent? If you answered yes, prove your claim. If you answered no, sketch an example which illustrates your claim, and outline a method that has the Newton method’s local convergence properties but which is globally convergent.

Numerical Methods for ODEs Consider the Linear Multistep Method \[ y_{n+2}-\frac {4}{3}y_{n+1}+\frac {1}{3}y_n=\frac {2}{3}hf_{n+2} \] for solving an initial value problem \(y'=f(y,x)\), \(y(0)=\eta \). You may assume that \(f\) is Lipschitz continuous with respect to \(y\) for all \(x\).

a) Analyze the consistency, stability, accuracy, and convergence properties of this method.

b) Would it be more reasonable to use this method or Euler’s method for the initial value problem: \[ y'=-10^8[y-\cos (x)]-\sin (x),\qquad y(0)=1? \] Justify your answer.

Discretization and Convergence. What is a discretization? What does it mean for a discretization to be convergent? Give a precise definition for a differential equation problem and a discretization of your choice.

Local Truncation Error. What is the local truncation error of a numerical scheme? What does consistency mean? Give a specific and precise example.

Stability. What does stability of a numerical scheme mean? Give a specific example.

Convergence is equivalent to stability and convergence. Describe and prove a specific example for this fundamental result.

Iterative Methods for Linear Systems Consider the boundary value problem \[ -u''(x)=f(x),\qquad \text {for }0\le x\le 1 \] with \(u(0)=u(1)=0\), and the following discretization of it: \[ -U_{j-1}+2U_j-U_{j+1}=F_j \] for \(j=1,2,\ldots ,N-1\) where \(Nh=1\) and \(F_j=h^2f(jh)\). Imagine applying the Jacobi iterative method to solve this system of linear equations. Analyze the convergence properties of the Jacobi method for this problem. In particular, express the speed of convergence as a function of the discretization stepsize \(h\). You may find it helpful to recall that the \((N-1)\times (N-1)\) matrix \(M\) defined by \[ m_{ij}= \begin {cases} 2, & \text {if } i=j,\\ -1, & \text {if } \lvert i-j\rvert =1,\\ 0, & \text {otherwise} \end {cases} \] has eigenvalues \[ \mu _l=4\sin ^2\left (\frac {\pi l}{2N}\right ),\qquad l=1,2,\ldots ,N-1. \] (NOTE: This is a model problem designed to give you a simple linear system on which to analyze the Jacobi method; one would certainly not use an iterative method to solve a tri-diagonal system in practice!)

Newton Type Methods Suppose \(f(x)\) is a real-valued function defined for all values of the real variable \(x\), and that \(f'\) is Lipschitz continuous. Consider the discrete Newton method: Guess \(x_0\) and for \(k=0,1,2,\ldots \), do: \[ a_k=\frac {f(x_k+h_k)-f(x_k)}{h_k} \] \[ x_{k+1}=x_k-\frac {f(x_k)}{a_k} \] Analyze the local convergence properties of this method. How does the speed of convergence depend on the choice of \(h_k\) (ignore any effect of round-off errors in this analysis)? Now taking roundoff into account what would be a reasonable size for \(h_k\) in practice?

Numerical Methods for ODEs Consider the Linear Multistep Method \[ y_{n+2}-\frac {4}{3}y_{n+1}+\frac {1}{3}y_n=\frac {2}{3}hf_{n+2} \] for solving an initial value problem \(y'=f(y,x)\), \(y(0)=\eta \). You may assume that \(f\) is Lipschitz continuous with respect to \(y\) for all \(x\). a) Analyze the consistency, stability, accuracy, and convergence properties of this method. b) Would this method be a reasonable one to use for the initial value problem: \[ y'=-10^8[y-\cos (x)]-\sin (x),\qquad y(0)=1? \] Justify your answer.

a) Derive a multistep method to approximate \(y'=f(t,y)\) of the form \[ y_{n+1}=y_n+h(b_{-1}y'_{n+1}+b_0y'_n+b_1y'_{n-1}) \] with \(y'_j=f(t_j,y_j)\), \(t_j=t_0+jh\) by interpolating \(f(s,y(s))\) at \(t_{n-1}\), \(t_n\), \(t_{n+1}\).

b) Analyze the local truncation error of this method and the consequences regarding convergence. What would be an appropriate predictor for this method?

c) Discuss the region of absolute stability of this method, and the predictor-corrector method described in part b).

Compare the stability of the Euler Predictor-Trapezoidal Corrector and Backwards Euler timestepping methods for solution of the heat equation \[ u_t=u_{xx}\qquad u(x,0)=u_0(x) \] using the centered second difference operator \[ u_j \to \frac {u_{j+1}-2u_j+u_{j-1}}{h^2}\qquad j=1,\ldots ,J \] where \(u_j\approx u(jh)\) on the interval \([0,2\pi ]\) with periodic boundary conditions, and mesh-size \(h=2\pi /J\). Discuss the accuracy of these methods?

a) Show that Jacobi’s method for solving \(Ax=b\) iteratively is convergent when \(A\) is strictly diagonally dominant. What does this imply about the spectral radius of another matrix obtained from \(A\).

b) Describe the successive over-relaxation (SOR) method with parameter \(\omega \) and it’s rate of convergence. Give a sufficient condition upon \(A\) and \(\omega \) for convergence.

c) Briefly compare the complexity of solving Poisson’s problem in three dimensions by direct and iterative methods. Assume the simplest methods, grid, differencing scheme, boundary conditions, etc.

a) Find the order/rate of convergence of Newton’s method to a multiple root \(r\) of multiplicity \(m>1\) of the scalar equation \(f(x)=0\).

b) What method is obtained by applying Aitken’s acceleration in this situation? Show that it is at least quadratically convergent.

c) How would you extrapolate an integration formula whose error behaved asymptotically as \[ I-I_n=c_3h^3+c_6h^6+\cdots \] where \(I\equiv \int _a^b f(x)\,dx\) is the exact value, \(I_n\) is the approximation with \(n\) subintervals, and \(h\equiv \frac {b-a}{n}\).