Project 1: Finite difference methods for the wave equation#

Topics covered: Finite differences, boundary conditions, stability, numerical error, convergence.

Report instructions#

  1. Your written report needs to include the derivations and plots that are requested under “Assignment 1-6” below. There are also a few questions that you need to answer.

  2. Submit your report as one pdf that includes the source code.

  3. In addition to the pdf, submit your source code files separately so that we can download them and test your code.

  4. All equations must be nicely formatted and all figures must have captions. The text must refer to the figures and further explain what they illustrate.

  5. The report does not need to have a specific structure (Introduction, Methods, Results, etc.). You are free to select your own headings as long as you include the requested derivations, plots, answers, and code.

  6. We recommend using LaTeX to generate the report. There are tutorials on e.g. Overleaf. It is also possible (but more cumbersome) to generate a nicely formatted report using e.g. Word.

  7. Put all group members’ names on the first page.

Introduction#

The acoustic wave equation can be written

\[ K^{-1} \phi_{tt} = \nabla \cdot \rho^{-1} \nabla \phi, \]

where \(K=K(\bar{x})\) is bulk modulus, \(\rho=\rho(\bar{x})\) is density, and \(\phi=\phi(\bar{x},t)\) is the (unknown) momentum potential, which satisfies

\[ p = - \phi_t, \quad \rho \bar{v} = \nabla \phi, \]

where \(p\) denotes pressure and \(\bar{v}\) the particle velocity vector. If density is constant, the acoustic wave equation can be written

\[ \phi_{tt} = c^2 \nabla^2 \phi, \quad ({\color{red}*}) \]

where \(c = \sqrt{K/\rho}\) is the speed of sound. The model \(({\color{red}*})\) is often called simply the wave equation and is often used to model different kinds of wave propagation (acoustic, electro-magnetic, etc.).

The 1D problem#

We start with the 1D case and restrict the computational domain to \(-1 \le x \le 1\). Since the system is on second order form (second derivative in time) two initial conditions are needed: \(\phi = \phi_0, \, \phi_t = \varphi_0\), \(-1 \le x \le 1\).

Consider the following IBVP,

\[\begin{split} \begin{array}{lll} \phi_{tt}=c^2\phi_{xx}, & -1 < x < 1 & t \in (0, T],\\ \alpha \phi_t + \phi_x=0, & x =-1, & t \in (0, T] ,\\ \beta \phi_t + \phi_x=0, & x =1, & t \in (0, T] ,\\ \phi = \phi_0, \quad \phi_t=\varphi_0, & -1 \le x \le 1, & t = 0,\\ \end{array} \end{split}\]

where \(c=c(x)>0\), and \(\alpha, \beta\), are some real constants.

Assignment 1#

Derive conditions on \(\alpha\) and \(\beta\) such that the IBVP is well posed. Show that the IBVP is energy-conserving in the special case of Neumann BC, given by \(\alpha=\beta=0\).

Spatial discretization in 1D#

The spatial domain (\(-1 \le x \le 1\)) is discretized using the following \(m+1\) equidistant grid points (with grid-spacing \(h\)):

\[ \begin{array}{lll} x_j=-1+jh,& j=0,1,\ldots ,m, & h=\frac{2}{m}\;. \end{array} \]

The approximate solution at grid point \(x_j\) is denoted \(\phi_j\), and the discrete solution vector is

\[ \boldsymbol{\phi} =[\phi_0, \phi_1, \ldots, \phi_m]^T. \]

A semi-discrete SBP approximation of the PDE (without BC or IC) can be written

\[\begin{split} \begin{array}{ll} \boldsymbol{\phi}_{tt}=c^2D_2\boldsymbol{\phi}, & t \ge 0,\\ \end{array} \end{split}\]

where \(D_2\) is a second-derivative SBP operator. After imposing homogeneous Neumann BC, the semi-discrete problem can be written in the form

\[ \boldsymbol{\phi}_{tt} = D \boldsymbol{\phi}, \]

for some matrix \(D\). To solve this second-order ODE using the fourth ourder Runge–Kutta method (RK4), note that one can re-write it as a first-order system by introducing the variable

\[\begin{split} \mathbf{w} = \begin{bmatrix} \boldsymbol{\phi} \\ \boldsymbol{\phi}_t \end{bmatrix}. \end{split}\]

Assignment 2#

Derive an SBP-SAT approximation of the IBVP and prove stability.

Exact solutions in 1D#

With Neumann BC and constant wave speed \(c\), the following standing wave is an exact solution to the IBVP (assuming that the initial data are chosen to match the exact solution!):

\[ \phi(x,t) = \cos(kx)\cos(\omega t), \]

where \(\omega = ck\) and \(k=n\pi\), \(n\in \mathbb{Z}\).

Assignment 3#

Write a Python or Matlab program that solves the SBP-SAT approximation derived in the previous assignment, for the special case of Neumann BC (\(\alpha=\beta=0\)). Use RK4 for time-integration.

Assignment 4#

Use your code from the previous assignment to perform a convergence study where you compare your discrete solution (\(\boldsymbol{\phi}\)) with the standing wave solution \(\phi(x,t)\) at the final time \(T=\pi\). Use the time-step restriction

\[ \Delta t \leq 0.1 \frac{h}{c}, \]

to ensure stability of the time-stepping. Note that you have to select the time step such that you do not overshoot the final time! We refer to the codes in Lab 1 for some inspiration. Use the parameter values \(c=3\) and \(k=2\pi\). Run the 2nd, 4th, and 6th order SBP operators for \(m=25,50,100,200,400\). Report the results in a “convergence plot”, i.e., a log-log plot of error as a function of \(h\). Use the discrete \(\ell^2\) norm of the error, which is defined as

\[ e_{\ell^2} = \sqrt{h} \sqrt{\sum \limits_{j=0}^m (\phi_j(T)-\phi(x_j,T))^2}. \]

Do you obtain the expected convergence behavior? Include appropriate reference lines corresponding to, for example, 2nd, 4th and 6th order convergence in your plot.

The 2D problem#

Now consider the wave equation in a rectangle of width 2 and height 1, with Neumann BC:

\[\begin{split}\label{eq:wave_PDE_2Da} \begin{array}{lll} \phi_{tt}= c^2 (\phi_{xx}+\phi_{yy}), & (x,y) \in \Omega, t>0\\ {\displaystyle \frac{\partial\phi}{\partial\hat{\mathbf{n}}} = 0}, & (x,y) \in \partial \Omega & t > 0,\\ \phi = f, & (x,y) \in \Omega, & t = 0,\\ \phi_t = g, & (x,y) \in \Omega, & t = 0,\\ \end{array} \end{split}\]

where \(\Omega = [-1, 1] \times [\frac{-1}{2}, \frac{1}{2}]\), \(\partial \Omega\) is the boundary of \(\Omega\), \(\hat{\mathbf{n}}\) is the outward unit normal, and \(\frac{\partial \phi}{\partial \hat{\mathbf{n}}}\) denotes the normal derivative of \(\phi\).

Assignment 5#

Use the energy method to show that the 2D IBVP is well posed.

Spatial discretization in 2D#

To simplify the 2D notation we will make use of the Kronecker product:

\[\begin{split} B \otimes C=\begin{bmatrix} b_{1,1}\,C &\cdots &b_{1,m}\,C\\ \vdots & & \vdots\\ b_{m,1}\,C &\cdots &b_{m,m}\,C\\ \end{bmatrix}\;, \end{split}\]

where \(B\) is an \(m \times m\) matrix and \(C\) is an \(n \times n\) matrix. The rectangular domain is discretized with an \(m_x \times m_y\)-point equidistant grid defined as

\[\begin{split} \begin{array}{lll} x_i=-1+ih_x,&i=0,1,\ldots,m_x,& h_x=\frac{2}{m_x},\\ y_j=-\frac{1}{2}+jh_y,& j=0,1,\ldots ,m_y,& h_y=\frac{1}{m_y}.\\ \end{array}\; \end{split}\]

The numerical approximation at grid point \((x_i,y_j)\) is denoted \(\phi_{i,j}\). The discrete solution vector is now defined as an \((m_x+1)(m_y+1) \times 1\) vector:

\[\begin{split} \begin{aligned} \boldsymbol{\phi}=[&\phi_{1,1},\, \phi_{2,1},\dots, \phi_{m_x,1},\\ &\phi_{1,2},\,\phi_{2,2},\dots , \phi_{m_x,2}, \\ &\vdots \\ &\phi_{1,m_y},\, \phi_{2,m_y},\dots \phi_{m_x,m_y}]^T. \end{aligned} \end{split}\]

Let \(D_m\) denote the discretisation matrix for the 1D problem with \(m+1\) grid points, with the SAT for Neumann BC included. Further, let \(I_n\) denote the \(n\times n\) identity matrix. The semi-discrete SBP-SAT approximation of the 2D problem can be written as

\[\begin{split} \begin{array}{ll} \boldsymbol{\phi_{tt}}=\left((I_{m_y+1} \otimes D_{m_x+1}) + (D_{m_y+1} \otimes I_{m_x+1})\right)\boldsymbol{\phi}, & t > 0. \quad ({\color{blue} *}) \\ \end{array} \end{split}\]

Assignment 6#

Write a Matlab or Python program that solves the 2D problem \(({\color{blue} *})\) using RK4 for time integration. Run a simulation with initial data

\[\begin{split} \begin{aligned} \phi(x,y,0) &= \exp\left(-\frac{(x-x_0)^2}{\sigma^2} - \frac{(y-y_0)^2}{\sigma^2} \right) \\ \phi_t(x,y,0) &= 0, \end{aligned} \end{split}\]

with \(x_0=y_0=0\) and \(\sigma = 0.05\). Further, set \(c=1\), \(m_x=200\), and \(m_y=100\), and use the 6th order SBP operators. Present a plot of the numerical solution at time \(T=3\).

Hint 1: To make the computations faster and lower the memory storage, make sure that you store the matrices in sparse format. We refer to the codes from Lab 1 for som inspiration.

Hint 2: To ensure that your 2D code is correct, you may want to test it on the following standing-wave solution:

\[ \phi(x,y,t) = \cos(k_x x)\cos(k_y y)\cos(\omega t), \]

where

\[ k_x = m\pi, \quad k_y = 2n\pi, \quad m,n \in \mathbb{Z}, \quad \omega = c \sqrt{k_x^2 + k_y^2}. \]

Use the standing-wave solution with, for example, \(k_x=k_y=2\pi\) to derive initial data for \(\phi\) and \(\phi_t\), and verify that the numerical solution matches the standing wave for \(t>0\). You do not need to report this test.