Math 5486: Final Exam
Time Allowed: 2 hours
Answer SECTION A (40%) and TWO of the three (30% + 30%) optional sections (B, C, D). If you have answered more than than the required two optional sections, you will be given credit for your best two optional sections. As a guide, I suggest you spend around 40 minutes on each section.
Please indicate your answers to multiple choice questions clearly with an ╳ to the left of your answers.
Please label any additional paper you use clearly with your name and page number.
Electronic devices are not needed and not permitted in this examination. You may use your own notes, limited to 2 sides of letter sized paper. Additional paper may be requested.
For grading purposes (leave blank)
A: Compulsory Section [40 points]
NoteUseful definitions
Consider the initial value problem: seek u:[0,T] \to \mathbb R such that
\begin{align} \tag{IVP} \begin{split} u(0) &= u_0 \\ u'(t) &= f\big(t,u(t)\big). \end{split} \end{align}
For a fixed mesh t_j = j \tau = \frac{jT}{n} (for 0 \leq j \leq n), we use the notation u_j \approx u(t_j) for approximations coming from the numerical scheme in question.
[10 points]
Please select all the true statements.
We can verify an exponential rate of convergence of an iterative scheme by plotting the iteration number against the error on a log-log plot.
All zeros of sparse matrices must be stored in memory and so it is difficult to work with large system sizes in practice.
There is f for which there are infinitely many u : [0,T] \to \mathbb R solving \text{(IVP)}.
One can obtain a third order accurate method for approximating \text{(IVP)} without knowledge of the derivatives of f.
Suppose a linear multi-step method has the characteristic polynomial \rho(z) = z(z-1). This method is zero-stable.
The region of absolute stability for backwards Euler is bounded.
The problem \text{(IVP)} with f(t,u) = u is well-posed.
For the Runge-Kutta method corresponding to the following Butcher tableau
\begin{align} \begin{array} {c|cccc} 0\\ c & a \\ \hline & b_1 & b_2 \end{array}\tag{A.1} \end{align}
to be consistent, we require b_1 + b_2 = 1.
Jacobi’s method for solving Ax = b is well-defined as long as A is symmetric and positive definite.
If \rho(T) < 1 (where \rho(T) is the spectral radius of T), then the iteration x^{(k+1)} = Tx^{(k)} + c converges to (I-T)^{-1}c.
In exact arithmetic the conjugate gradient method converges to the exact solution to Ax = b in at most n steps where A\in\mathbb R^{n\times n} is symmetric and positive definite.
2.
[5 points]
Show that Euler’s method is first order accurate.
3.
[5 points]
When using an adaptive Runge-Kutta method, one uses an estimate for the local error so that one may refine the grid where this error estimator is large. Why do we not just use the actual error instead?
4.
[10 points]
Suppose A \in \mathbb R^{n\times n} is symmetric and positive definite and consider the function g(x) := \frac{1}{2} x^T A x - b^T x. Show that, if x is a minimiser of g, then Ax = b.
Hint: you may use the fact that, if x is a minimiser of g, then
\begin{align} 0 &\leq \frac{g(x+\epsilon y)-g(x)}{\epsilon} % = y^T (A x - b) + \frac{\epsilon}{2} y^T A y. \tag{A.2} \end{align}
for all y\in \mathbb R^n and \epsilon> 0.
5.
[10 points]
Define the mesh x_i = a + i h where h := \frac{b-a}{m} and i=0,\dots,m. Suppose that u solves the boundary value problem
\begin{align} &\mathcal L[u] = f(x) \quad \text{on }(a,b) \nonumber\\ &u(a) = u(b) = 0 \tag{A.3} \end{align}
and U_h solves a finite difference approximation given by
\begin{align} &\mathcal L_h U_h = F \nonumber\\ &(U_h)_0 = (U_h)_m = 0 \tag{A.4} \end{align}
where F_j := f(x_j). Define (I_h u)_j := u(x_j) to be the restriction of u to the mesh. Suppose the finite difference approximation is consistent and stable: there exists p>0 and C_1, C_2 >0 such that
\begin{align} \| \mathcal L_h (I_h u) - F \| &\leq C_1 h^p, \tag{consistency}\\ % \| V - W \| &\leq C_2 \| \mathcal L_h V - \mathcal L_h W \| \tag{stability} \end{align}
for all V, W. Show that \| U_h - I_h u \| \leq C_1 C_2 h^p (i.e. the method is convergent of order p).
Choose TWO sections from B, C, D (each worth 30 points)
B: FEM: Poisson’s equation [30 points]
In this question, we shall consider Poisson’s equation on [a,b]. In class, we mainly focused on homogeneous Dirichlet and periodic boundary conditions. Here, we look at inhomogeneous Dirichlet and homogeneous Neumann boundary conditions. We consider the mesh x_i := a + i\frac{b-a}{n} =: a + i h for i=0,\dots,n and the corresponding hat functions \phi_0,\dots,\phi_n. Recall that \phi_i is linear on [x_j,x_{j+1}] for all j=0,\dots,n-1 and \phi_i(x_j) = \delta_{ij} (i.e. \phi_i(x_i) = 1 and \phi_i(x_j)=0 for all i\not=j). In this question, you may use without proof that the discrete Laplacian is given by
\begin{align} - L_{ij} &:= \int_a^b \phi_i'(x) \phi_j'(x) \mathrm{d}x % = \frac{1}{h} \begin{cases} 2 &\text{if } i=j \nonumber\\ -1 &\text{if } i=j\pm1 \nonumber\\ 0 &\text{otherwise} \end{cases} \tag{B.1} \end{align}
for i,j = 1,\dots,n-1.
We first consider Poisson’s equation with Dirichlet boundary conditions:
\begin{align} \begin{split} -&u''(x) = f(x) \quad \text{on }(a,b) \nonumber\\ &u(a) = A, \quad u(b) = B. \end{split} \tag{B.2} \end{align}
- Show that, if u solves \text{(B.2)}, then
\begin{align} \int_{a}^b u'(x) \phi'_i(x) \mathrm{d}x = \int_{a}^b f(x) \phi_i(x) \mathrm{d}x \tag{B.3} \end{align}
\quad\,\, for all i = 1,\dots,n-1. [5 points]
Now approximate the solution of \text{(B.2)} with u_h(x) := U_0 \phi_0(x) + \sum_{j=1}^{n-1} U_j \phi_j(x) + U_n \phi_n(x) for some coefficients (U_j).
What should you define U_0 and U_n as in order for the approximate solution u_h to satisfy the boundary conditions in \text{(B.2)}? [5 points]
Suppose that u_h satisfies \text{(B.3)} and write down the linear system of equations that is satisfied by (U_j)_{j=1}^{n-1}. [5 points]
Now, consider Poisson’s equation with homogeneous Neumann boundary conditions:
\begin{align} \begin{split} -&v''(x) = f(x) \quad \text{on } (a,b) \nonumber\\ &v'(a) = 0, \quad v'(b) = 0. \end{split} \tag{B.4} \end{align}
with \int_a^b f(x)\mathrm{d}x = 0.
- Show that, if v is a solution to \text{(B.4)}, then v + C is also a solution for all constants C. [5 points]
If v solves \text{(B.4)}, then, by integrating by parts, we obtain
\begin{align} \int_a^b v'(x) \phi_i'(x) \mathrm{d}x = \int_a^b f(x) \phi_i(x) \mathrm{d}x \tag{B.5} \end{align}
for all i=0,\dots,n. We approximate v with v_h(x) := \sum_{j=0}^{n} V_j \phi_j(x) for some coefficients (V_j)_{j=0}^{n} (notice we now include the “half-hat” functions \phi_0 and \phi_n).
- Write down the linear system of equations -L_{\rm N} V = F satisfied by (V_{j})_{j=0}^n if v_h satisfies \text{(B.5)}. Here, L_{\rm N} is the discrete Neumann Laplacian. [5 points]
It turns out that -L_{\rm N} is singular which reflects the fact that the exact solution is only defined up to a constant C. In practice, we “pin” the initial point V_0 = C by replacing the first row of -L_{\rm N} with \begin{pmatrix}1, 0,\dots, 0\end{pmatrix} and the first entry of F with F_0 := C. Then, solving the resulting linear system of equations gives an approximation to the exact solution with v(0) = C. Here, we plot the FEM solutions for both Dirichlet (A = \frac{1}{50}, B = -\frac{1}{50}) and homogeneous Neumann boundary conditions:
- Which one of these numerical solutions approximates the Neumann problem \text{(B.4)}? (You can just answer “top” or “bottom”). [2 points]
- Is the convergence exponential, O(e^{-c\,n}) for c > 0, or algebraic, O(n^{-p}) for p > 0? [3 points]
Choose TWO sections from B, C, D (each worth 30 points)
C: Method of Lines: Transport equation [30 points]
Suppose we are interested in solving the transport equation
\begin{align} \begin{split}&u_t(x,t) = u_x(x,t) \text{ for } x\in(a,b), t>0 \nonumber\\ &u(x,0) = u_0(x) \nonumber\\ &u(a,t) = u(b,t) \end{split}\tag{C.1} \end{align}
for some initial condition u_0(x).
We again take x_i := a + \frac{b-a}{m} i =: a + i h (for i=0,\dots,m) but now approximate the solution at time t with the (time-dependent) vector U(t):
\begin{align} U(t) = \begin{pmatrix} U_1(t) \\ U_2(t) \\ \vdots \\ U_{m}(t) \end{pmatrix} % \approx \begin{pmatrix} u(x_1,t) \\ u(x_2,t) \\ \vdots \\ u(x_{m}, t) \end{pmatrix}. \tag{C.2} \end{align}
This is known as the semi-discretisation step because we have discretised in space but not time. Notice, we are considering periodic boundary conditions so I have omitted U_{0}(t) = U_m(t) from this vector.
- Show that approximating u(x_i,t) with the centered difference \frac{U_{i+1}(t) - U_{i-1}(t)}{2h} results in the following system of differential equations
\begin{align} U'(t) &= D^0 U(t) \nonumber\\ U(0) &= \big( u_0(x_i) \big)_{i=1}^m \tag{C.3} \end{align}
\quad\,\, where D^0 \in \mathbb R^{m\times m} is given by
\begin{align} D^0 &= \frac{1}{2h} \begin{pmatrix} 0 & 1 & & -1\\ -1 & 0& \ddots & \\ & \ddots& \ddots & 1 \\ 1 & & -1 & 0 \end{pmatrix} \tag{C.4} \end{align}
\quad\,\, (where the upper diagonal entries are all 1 and the lower diagonal entries are all -1). [5 points]
Since D^{0} is a normal matrix, it has an eigenvalue decomposition D^0 = P \Lambda P^{-1} (where \Lambda is diagonal and consists of eigenvalues of D^0 and the columns of P are the corresponding eigenvectors).
- Show that, if U(t) solves \text{(C.3)}, then W(t) := P^{-1}U(t) solves the (decoupled) equations:
\begin{align} W'(t) = \Lambda W(t). \tag{C.5} \end{align}
\quad\,\, [5 points]
Thus, we have W_i'(t) = \lambda_{i} W_i(t) for i=1,\dots,m where \lambda_i := \Lambda_{ii} are the eigenvalues of D^0.
- Show that the exact solution to W_i'(t) = \lambda_i W_i(t) can be written W_i(t) = W_i(0) \mathrm{e}^{\lambda_i t}. [5 points]
Therefore, if \mathrm{Re} \, \lambda_i \leq 0, then the solution to \text{(C.5)} and thus the solution to \text{(C.3)} remains bounded as t\to \infty.
Define the Fourier basis: (\bm e_k)_j := \frac{1}{\sqrt{m}} \mathrm{e}^{\frac{2\pi \mathrm{i}jk}{m} } for k=1,\dots,m.
- Show that D^0 \bm e_k = \frac{\mathrm{i}}{h} \sin(\frac{2\pi k}{m}) \bm e_{k}. [5 points]
You have therefore seen that the eigenvalues, \lambda_k, of D^0 are purely imaginary and |\lambda_k| \leq \frac{1}{h} and W(t) remains bounded as t\to\infty. Now consider a grid t_j := j \tau and approximate W_i(t_j) \approx W_{i,j} with Euler’s method:
\begin{align} W_{i,j+1} := W_{i,j} + \tau \lambda_i W_{i,j}. \tag{C.6} \end{align}
- Explain why this leads to an unstable numerical method. [5 points]
Instead one can consider Backward Euler:
\begin{align} W_{i,j+1} := W_{i,j} + \tau \lambda_i W_{i,j+1}. \tag{C.7} \end{align}
- Explain why this leads to a stable method. [5 points]
TURN OVER FOR SECTION D
Choose TWO sections from B, C, D (each worth 30 points)
D: Finite differences: Heat equation [30 points]
Let us consider the following heat equation
\begin{align} \begin{split} &u_t(x,t) = u_{xx}(x,t) \nonumber\\ &u(x,0) = u_0(x) \nonumber\\ &u(a,t) = u(b,t) = 0. \end{split}\tag{D.1} \end{align}
For 0 \leq \theta \leq 1, we consider the following finite difference approximation to the heat equation:
\begin{align} \frac{u_{i,j+1}-u_{i,j}}{\tau} = % (1-\theta) \frac{u_{i-1,j}-2u_{i,j}+u_{i+1,j}}{h^2} % + \theta \frac{u_{i-1,j+1}-2u_{i,j+1}+u_{i+1,j+1}}{h^2} \tag{D.2} \end{align}
where x_{i} = a + i h and t_j = j \tau. Recall, we hope that u_{i,j} will approximate the exact solution to \text{(D.1)} at (x_i, t_j).
- For which \theta do we recover the \text{(FTCS)} (forward in time, centered in space) and \text{(BTCS)} (backward in time, centered in space) methods that we saw in lectures? [5 points]
Recall that, both \text{(FTCS)} and \text{(BTCS)} are consistent of order O(\tau + h^2).
- Use the Taylor expansions
\begin{align} \frac{u(x,t+\tau)-u(x,t)}{\tau} % &= u_t(x,t) + \frac{\tau}{2} u_{tt}(x,t) + O(\tau^2) \nonumber\\ % \frac{u(x-h,t)-2u(x,t)+u(x+h,t)}{h^2} % &= u_{xx}(x,t) + \frac{h^2}{12}u_{xxxx}(x,t) + O(h^4) \tag{D.3} \end{align}
\quad\,\, to show that \text{(D.2)} is consistent. [5 points]
- For what \theta does \text{(D.2)} have order of accuracy O(\tau^2 + h^2). Hint: You will need to use u_{xxt} = (u_{xx})_t = (u_{t})_t = u_{tt} if u solves the equation. [10 points]
Recall that \text{(FTCS)} is stable when \mu := \frac{\tau}{h^2} \leq \frac{1}{2} and \text{(BTCS)} is unconditionally stable.
- Show that \text{(D.2)} is von Neumann stable for \theta \geq \frac{1}{2}. Hint: consider u_{i,j} := \mathrm{e}^{\mathrm{i} k x_i} and u_{i,j+1} := \lambda(k) \mathrm{e}^{\mathrm{i} k x_i} and use the identity 1-\cos x = 2\sin^2 \frac{x}2. [10 points]
END