Finite Difference Methods

(Finite difference formulas) Given \(h > 0\), compute weights \(\{w_j\}_{j=0}^2\) for the following one-sided finite difference formula for the first derivative \(f'(x)\): \[ f'(x) \approx \sum _{j=0}^2 w_j f(x+jh) = w_0 f(x) + w_1 f(x+h) + w_2 f(x+2h). \]

What is the order of accuracy of your formula?

Finite difference formulas

Given \(h>0\), compute weights \(v_0\) and \(\{w_j\}_{j=-1}^{1}\) for the following finite difference formula for the third derivative \(f^{(3)}(x)\): \[ f^{(3)}(x)\approx v_0f'(x)+\sum _{j=-1}^{1}w_j f(x+jh) \] \[ =v_0f'(x)+w_{-1}f(x-h)+w_0f(x)+w_1f(x+h). \] What is the order of accuracy of your formula?

Higher Order Methods

Consider the following problem for \(v(x)\) on \([0,1]\): \[ v''(x)=f(x), \tag {2} \] with \(v(0)=v(1)=0\). Let \(N\cdot h=1\) and define \(x_j=j\cdot h\) for \(j=0,1,\ldots ,N\). The finite-difference scheme \[ \Delta _h U_j^h\equiv \frac {1}{h^2}\left (U_{j-1}^h-2U_j^h+U_{j+1}^h\right ) = f_j^h, \] for \(j=1,\ldots ,N-1\) with \(U_0^h=U_N^h=0\) gives values \(U_j^h\) that approximate \(v(x_j)\) with an error of \(O(h^2)\). Here the superscript \(h\) is used to indicate the grid size for the solution. Show how to use this method to find a numerical solution \(W_j\) whose values approximate \(v(x_j)\) with an error of \(O(h^4)\) and a numerical solution \(Y_j\) whose values approximate \(v(x_j)\) with an error of \(O(h^6)\).

Finite difference formulas

Given \(h>0\), compute weights \(w_j\) for the following 5-point finite difference formula for the third derivative \(f^{(3)}(x)\): \[ f^{(3)}(x)\approx \sum _{j=-2}^{2}w_j f(x+jh) \] \[ = w_{-2}f(x-2h)+w_{-1}f(x-h)+w_0f(x)+w_1f(x+h)+w_2f(x+2h). \] What is the order of accuracy of your formula?

Numerical ODEs

a) The difference equation \[ y_{n+k}+a_1y_{n+k-1}+\cdots +a_ky_n=0 \] can, after the substitution \(\mathbf {y}_j=(y_j,y_{j+1},\ldots ,y_{j+k-1})^T\), be written in the form \(\mathbf {y}_{n+1}=A\mathbf {y}_n\). What does the matrix \(A\) look like? Describe the condition for \(y\) to remain bounded in terms of \(A\).

b) The differential equation \[ y^{(k)}+a_1y^{(k-1)}+\cdots +a_ky=0 \] can, after the substitution \(\mathbf {y}=(y,y',\ldots ,y^{(k-1)})^T\), be written in the form \(\mathbf {y}'=A\mathbf {y}\). What does the matrix \(A\) look like? Describe the condition for \(y\) to remain bounded in terms of \(A\).

c) Show that all solutions of the difference equation \[ y_{n+1}-2\lambda y_n+y_{n-1}=0 \] are bounded, as \(n\to \infty \), if \(\lambda \) is real and \(\lvert \lambda \rvert <1\).

Finite-Difference Methods Consider the initial boundary value problem \[ v_t=\beta v_{xx},\qquad 0<x<1,\quad t>0,\qquad v(x,0)=f(x) \] with Dirichlet boundary conditions \[ v(0,t)=0,\qquad v(1,t)=0. \] Suppose that the interval \([0,1]\) is partitioned by meshpoints \(x_j=jh\) for \(j=0,1,\ldots ,N\) and \(h=1/N\). A possible discretization of the above problem uses the Crank-Nicolson scheme for the differential equation: \[ \left (I-\frac {\beta k}{2}D_+D_-\right )u_j^{n+1} = \left (I+\frac {\beta k}{2}D_+D_-\right )u_j^n \] for \(j=1,\ldots ,N-1\), combined with the boundary conditions: \[ u_0^n=0\qquad \text {and}\qquad u_N^n=0 \] for each time level \(n\). Initial data is obtained by setting \(u_j^0=f_j\) for \(j=0,1,\ldots ,N\). (Note that \(D_+\) and \(D_-\) are the usual forward and backward difference operators.) Set up the matrix system that must be solved in order to advance the numerical solution from time level \(n\) to time level \(n+1\). Show that the suggested discretization is stable in the \(\ell _2\) norm \[ \lVert u\rVert _2=\left (\sum _{j=0}^N h\lvert u_j\rvert ^2\right )^{1/2} \] for all choices of timestep \(k\) and spacestep \(h\). (Note that von Neumann analysis is not enough to prove stability for the initial boundary value problem.) What other issues are important to consider in deciding whether this scheme is a good choice for the problem?

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_{ij}= \begin {cases} -2, & i=j,\\ 1, & \lvert i-j\rvert =1,\\ 0, & \text {otherwise} \end {cases} \qquad (1\le i,j\le N-1) \] 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 \(M\) defined by \[ M_{ij}= \begin {cases} -1, & i=j\in \{0,N\},\\ -2, & i=j,\ 1\le i\le N-1,\\ 1, & \lvert i-j\rvert =1,\\ 0, & \text {otherwise} \end {cases} \qquad (0\le i,j\le N) \] 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\).

Iterative Methods for Linear Systems

a) (6%) Consider the linear system \(Ax=b\) and a splitting of \(A:A=M-K\), with \(M\) being nonsingular, an iterative method is to iterate \[ 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.

b) (14%) To solve the boundary value problem \[ -u''(x)=f(x),\qquad \text {for }0<x<1, \] with \(u'(0)=0\) and \(u(1)=0\), we use second-order central differences for the discretization of the equation at mesh points \(x_j=j/N\): \[ -u_{j-1}+2u_j-u_{j+1}=f_j=\frac {f(x_j)}{N^2}, \qquad j=0,1,\ldots ,N-1 \] and the central difference approximation for \(u'\) at \(x_0=0\): \(u_1-u_{-1}=0\), to eliminate the ghost point.

Show that the finite difference equations and the boundary conditions lead to the following linear system \(Au=b\), where \[ \setlength {\arraycolsep }{5pt} A= \begin {bmatrix} 1 & -1 & 0 & \cdots & 0 & 0\\ -1 & 2 & -1 & \cdots & 0 & 0\\ 0 & -1 & 2 & \ddots & 0 & 0\\ \vdots & \vdots & \ddots & \ddots & -1 & 0\\ 0 & 0 & 0 & -1 & 2 & -1\\ 0 & 0 & 0 & 0 & -1 & 2 \end {bmatrix}, \quad u= \begin {pmatrix} u_0\\ u_1\\ \vdots \\ u_{N-2}\\ u_{N-1} \end {pmatrix}, \quad b=\frac {1}{N^2} \begin {pmatrix} \frac {1}{2}f_0\\ f_1\\ \vdots \\ f_{N-2}\\ f_{N-1} \end {pmatrix}. \] Notice that equations are scaled so that \(A\) is symmetric. Determine if \(A\) has eigenvectors \[ v^{(k)}= \begin {bmatrix} 1\\ \cos \frac {1}{2}k\pi x_1\\ \vdots \\ \cos \frac {1}{2}k\pi x_j\\ \vdots \\ \cos \frac {1}{2}k\pi x_{N-2}\\ \cos \frac {1}{2}k\pi x_{N-1} \end {bmatrix}, \qquad k=1,3,5,\ldots ,2N-3,2N-1, \] and use them to find the eigenvalues of the iteration matrix \(R\) for the Jacobi method. Analyze the convergence properties of the Jacobi method applied to this problem, and 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\)? In reality, should you use an iterative method to solve this matrix problem? If not, what other method is preferred?

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.

Iterative Methods for Linear Systems

(a) Consider the linear system \(Ax=b\) and a splitting of \(A:A=M-K\), with \(M\) being nonsingular, an iterative method is to iterate \[ 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.

(b) To solve the boundary value problem \[ -u''(x)=f(x),\qquad \text {for }0<x<1, \] with \(u(0)=0\), and \(u'(1)=0\), we use second-order central differences for the discretization at the mesh points \(x_j=j/N\): \[ -u_{j-1}+2u_j-u_{j+1}=f_j=\frac {f(x_j)}{N^2}, \qquad j=1,2,\ldots ,N \] and the central difference for \(u'\) at \(x_N\) to eliminate the ghost point \(x_{N+1}\): \[ u_{N+1}=u_{N-1}. \] Use all the finite difference equations and the boundary conditions to obtain the linear system \(Au=b\), where \[ \setlength {\arraycolsep }{5pt} A= \begin {bmatrix} 2 & -1 & 0 & \cdots & 0 & 0\\ -1 & 2 & -1 & \cdots & 0 & 0\\ 0 & -1 & 2 & \ddots & 0 & 0\\ \vdots & \vdots & \ddots & \ddots & -1 & 0\\ 0 & 0 & 0 & -1 & 2 & -1\\ 0 & 0 & 0 & 0 & -2 & 2 \end {bmatrix}, \quad u= \begin {pmatrix} u_1\\ u_2\\ \vdots \\ u_{N-1}\\ u_N \end {pmatrix}, \quad b=\frac {1}{N^2} \begin {pmatrix} f_1\\ f_2\\ \vdots \\ f_{N-1}\\ f_N \end {pmatrix}. \]

Verify that \(A\) has eigenvectors \[ v^{(k)}= \begin {bmatrix} \sin \frac {k\pi }{2}x_1\\ \sin \frac {k\pi }{2}x_2\\ \cdot \\ \cdot \\ \sin \frac {k\pi }{2}x_j\\ \cdot \\ \cdot \\ \sin \frac {k\pi }{2}x_{N-1}\\ \sin \frac {k\pi }{2}x_N \end {bmatrix}, \qquad k=1,3,5,\ldots ,2N-3,2N-1, \] and use them to find the eigenvalues of the iteration matrix \(R\) for the Jacobi method. Analyze the convergence properties of the Jacobi method applied to this problem, and 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\)?

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).

Finite-Difference Methods Consider the initial boundary value problem \[ v_t=\beta v_{xx},\qquad 0<x<1,\qquad 0<t,\qquad v(x,0)=f(x) \] with Dirichlet boundary conditions \[ v(0,t)=0,\qquad v(1,t)=0. \]

Suppose that the interval \([0,1]\) is partitioned by meshpoints \(x_j=jh\) for \(j=0,1,\ldots ,N\) and \(h=1/N\). A possible discretization of the above problem uses the Crank-Nicolson scheme for the differential equation: \[ \left (I-\frac {\beta k}{2}D_+D_-\right )u_j^{n+1} = \left (I+\frac {\beta k}{2}D_+D_-\right )u_j^n \] for \(j=1,\ldots ,N-1\), combined with the boundary conditions: \[ u_0^n=0\qquad \text {and}\qquad u_N^n=0 \] for each time level \(n\). Initial data is obtained from by setting \(u_j^0=f_j\) for \(j=0,1,\ldots ,N\). (Note that \(D_+\) and \(D_-\) are the usual forward and backward difference operators.)

Set up the matrix system that must be solved in order to advance the numerical solution from time level \(n\) to time level \(n+1\). Show that the suggested discretization is stable in the \(l_2\) norm \[ \lVert u\rVert _2= \left (\sum _{j=0}^N h\lvert u_j\rvert ^2\right )^{1/2} \] for all choices of timestep \(k\) and spacestep \(h\). What other issues should be considered in deciding whether the suggested scheme is a good choice for the problem?

Finite-Difference Methods Consider the spatially-periodic initial value problem for \(v(x,t)\): \[ v_t=cv_x,\qquad \text {for }0\le x\le 1 \] with \[ v(x,0)=f(x),\qquad \text {for }0\le x\le 1. \] Here the function \(f\) has period 1, and \(c\) is a positive constant. Note that the exact solution of the problem is \(v(x,t)=f(x+ct)\).

Consider the following difference approximations to the PDE: \[ \begin {aligned} \text {i)}\quad &u_j^{n+1}=u_j^n+\frac {ck}{h}(u_{j+1}^n-u_j^n),\\ \text {ii)}\quad &u_j^{n+1}=u_j^n+\frac {ck}{h}(u_j^n-u_{j-1}^n),\\ \text {iii)}\quad &u_j^{n+1}=u_j^n+\frac {ck}{2h}(u_{j+1}^n-u_{j-1}^n). \end {aligned} \] In each case, \(Nh=1\), and the difference scheme is supplemented with the discrete initial conditions \(u_j^0=f(jh)\) for \(j=0,1,\ldots ,N\) and boundary conditions \(u_N^n=u_0^n\), \(u_{-1}^n=u_{N-1}^n\), and \(u_{N+1}^n=u_1^n\) for all time levels \(n\ge 0\). For each scheme, tell whether it is a reasonable scheme to use on the above problem, and justify your answers.