Lab 1: Finite difference methods#

Topics covered: Well-posedness, stability, numerical error, convergence.

Note: The labs are not mandatory, but highly recommended.

Getting started#

Download the following Python files (right-click the links and click Download linked file):

Lab11.py, operators.py, rungekutta4.py

Open the the file Lab11.py and try to run it. You should see a figure window with a moving pulse. If the code runs, you can move on to the next section.

Note: The code uses the modules numpy, scipy and matplotlib, which may not be installed on your computer. On Mac/Unix, you can install e.g. numpy by typing the following in a terminal window:

python3 -m pip install numpy

The University’s Windows computers have VS Code installed. Open VS Code and press ctrl+shift+P. Search for python interpreter and click Python: Select Interpreter in the drop-down menu that appears. Make sure that the newest version (probably 3.10.x) is selected. Then, in the terminal window, type the following at the prompt:

& "C:/Program Files/Python310/python.exe" -m pip install numpy.

Repeat for the other modules, as required.

Part 1: the advection equation#

We start by considering the advection equation with periodic boundary conditions:

\[\begin{split}\begin{align} u_t + cu_x &= 0, \quad x \in [0, L] \\ u(0,t) &= u(L, t) \\ u(x, 0) &= f(x) \end{align} \end{split}\]

For simplicity we will consider the domain length \(L=1\) and the wave speed \(c=1\), in appropriate units. We will also assume that the initial data function \(f\) is periodic with period \(L\). The exact solution to the above problem is

\[ u(x,t) = f(x-ct) . \]

In Lab11.py we consider the Gaussian-bell initial data:

\[ f(x) = e^{-{\left( \frac{x-0.15}{0.05} \right)^2}}. \]

The code uses standard centered finite difference approximations of orders 2, 4, 6, 8, 10, and 12 to solve the PDE numerically. The grid spacing is \(h = 1/m\), where \(m\) is the number of grid intervals. You can change the order of accuracy and/or \(m\) in the main() method in Lab11.py.

The code prints the \(l^2\) error at the final time \(t=3\), compared to the exact solution. Let

\[ \mathbf{v} = [v_0, v_1, \ldots, ..., v_{m-1}] \]

denote the discrete solution vector, where \(v_j\) is the discrete approximation at gridpoint \(x_j\). The \(l^2\) error \(e\) is defined as

\[ e = \sqrt{h\sum_j (v_j-u(x_j))^2}. \]

Dispersion errors#

  1. Run Lab11.py for some different values of \(m\) (e.g. 50, 100 and 200) and different orders of accuracy. Does the approximation seem to improve if \(m\) is increased?

  2. Set \(m=50\) and try all the different orders. Does the approximation seem to improve with increasing order?

The advection equation translates all wavenumbers with the same speed (hence the initial pulse retains its shape). However, the numerical approximation gets the wave speed wrong, in particular for large wavenumbers. The error casued by the difference in wave speed is known as dispersion error. For long-time simulations even a small difference in wave speed will lead to a large error at the final time.

Convergence#

As mentioned, the code uses finite differences to approximate the spatial derivatve. For time integration, the code uses the classical 4th order Runge–Kutta method, with the time step selected as

\[ \Delta t = C \frac{h}{c}. \]

The non-dimensional constant \(C\) is here chosen as \(C=0.1\)

  1. For orders 2, 4, and 6, run the code with \(m=100\), \(m=200\), and \(m=400\). How much (by what factor) does the error decrease when \(m\) is doubled? Is this the convergence rate that you would expect from theory? If not, why might that be? Bonus task: Write a script that runs these simulations and prints the errors and convergence rates.

  2. Try the same \(m\) values with the 12th order finite difference method. Try to explain the convergence rates that you observe!

Programming#

Go through the code and try to understand how it works. The codes provided in the labs may serve as a helpful starting point when solving the project assignments.

Part 2: Maxwell’s equations#

For this part you will additionally need the following programs: Maxwell_1D_BC.py, Maxwell_1D_Interface.py, Maxwell_2D_Interface.py

Assuming no free charges or currents, the 2D Maxwell’s equations can be written (after appropriate dimensional scalings)

\[ {\bf C} {\bf u}_t={\bf A} {\bf u}_x+{\bf B} {\bf u}_y\;, \]

where

\[\begin{split} {\bf u}= \begin{bmatrix} E^{(x)}\\ H^{(z)}\\ E^{(y)}\\ \end{bmatrix} \end{split}\]

denotes the components of the electric and magnetic fields. (Here \(E^{(x)}\) is for example the electric field in the x-direction and \(H^{(z)}\) the magnetic field in the z-direction. ) The matrices \({\bf A}\), \({\bf B}\) and \({\bf C}\) are given by

\[\begin{split} {\bf A}= \begin{bmatrix} 0& 0& 0\\ 0& 0& -1\\ 0& -1& 0\\ \end{bmatrix} \,,\, {\bf B}= \begin{bmatrix} 0& 1& 0\\ 1& 0& 0\\ 0& 0& 0\\ \end{bmatrix} \,,\, {\bf C}= \begin{bmatrix} \epsilon& 0& 0\\ 0& \mu& 0\\ 0& 0& \epsilon\\ \end{bmatrix}\;. \end{split}\]

The permittivity \(\epsilon\) and permeability \(\mu\) may vary in space, but are here assumed to be time-independent. Many of the difficulties that arise in finite difference approximation of Maxwells’s equations are related to material interfaces and boundary conditions.

Boundary conditions in 1D#

We first consider Maxwell’s equations in 1D. Removing the \(E^{(x)}\) component leads to the model

\[ {\bf C} {\bf u}_t={\bf A} {\bf u}_x\;, \]

where now

\[\begin{split} {\bf u}= \begin{bmatrix} E\\ H\\ \end{bmatrix} \,,\, {\bf A}= \begin{bmatrix} 0& 1\\ 1& 0\\ \end{bmatrix} \,,\, {\bf C}= \begin{bmatrix} \epsilon& 0\\ 0& \mu\\ \end{bmatrix}\;. \end{split}\]

\(E\) and \(H\) denote the (mutually perpendicular) tangential electric and magnetic field components.

There are many possible boundary conditions for Maxwell’s equations. Physics usually tells us which conditions to impose, but it is important to check that the resulting IBVP is well-posed. Loosely speaking, well-posedness requires that we impose the minimal number of boundary conditions that leads to a bounded solution. Here, we will test 4 different sets of boundary conditions:

  1. Specifying \(E\) at the two boundaries.

  2. Specifying \(H\) at the two boundaries.

  3. Specifying a linear combination \(E+\beta\,H\) and \(E-\beta\,H\) at the left and right boundary, respectively.

  4. Specifying both \(E\) and \(H\) at the two boundaries.

Run Maxwell_1D_BC.py. The numerical solution is integrated in time to \(t=3.6\) using the fourth order accurate Runge-Kutta method and time-step \(\Delta t =\frac{h}{10}\), to make the time-integration errors negligible. You can specify the number of grid points \(m\) (i.e, the resolution in space) and the order of accuracy of the SBP finite difference operator in the code. Try the 4 types of boundary conditions listed above by changing the parameter bc_type in the code. For boundary condition type 3, you can also specify the value of \(\beta\). Try different values of \(\beta\) including \(\pm 1\) and \(\pm 2\). In the special case \(\beta=0\) this condition reduces to condition 1 (specifying the electric component). In the numerical tests it is recommended to use at least \(m=101\) grid-points.

Questions:

  1. What happens when \(\beta\) is positive? Hint: BC set 3 is ill-posed for \(\beta>0\).

  2. What happens for \(\beta=-1\)? For \(\beta=-2\)? Try to explain. Hint: The characteristic variables of this system are \(E \pm \sqrt{\mu/\epsilon} H\).

  3. What happens when you specify boundary condition type 4, using for example SBP operators of order \(6\)? Hint: BC set 4 is ill-posed.

  4. How does ill-posedness manifest in case 3 compared to case 4? Try to explain the difference.

Coupling of two different media in 1D#

Consider Maxwell’s equations in 1D in the presence of a dielectric boundary, i.e., an internal boundary between two different materials. The resulting system can be written,

\[\begin{split} \begin{array}{ll} {\bf C}^{(l)}{\bf u}^{(l)}_t={\bf A} {\bf u}^{(l)}_x, & x \in [x_l,\,x_I] \\ {\bf C}^{(r)}{\bf u}^{(r)}_t={\bf A} {\bf u}^{(r)}_x, & x \in [x_I,\,x_r] \\ \end{array} \end{split}\]

where

\[\begin{split} {\bf C}^{(l,\, r)}= \begin{bmatrix} \epsilon^{(l,\, r)}& 0\\ 0& \mu\\ \end{bmatrix}\;, \end{split}\]

and \(\epsilon^{(l)} \ne \epsilon^{(r)}\). Hence, at the media interface \(x=x_I\) there is a discontinuity in permittivity. The wave speeds of the two media are

\[ c^{(l,r)}=\frac{1}{\sqrt{\epsilon^{(l,r)}\mu}} \]

and the wave impedances are

\[ Z^{l,r} = \sqrt{\frac{\mu}{\epsilon^{(l,r)}}}, \]

Assume that we use a Gaussian profile in the left sub-domain as initial data for the electric and magnetic fields. The Gaussian profile will split into left- and right-going Gaussian waves propagating with speed \(c^{(l)}\). The right-going wave eventually hits the media interface at \(x=x_I\) and is further divided into a reflected wave and a transmitted wave. The transmission (\(T\)) and reflection (\(R\)) coefficients for the electric component are given by

\[\begin{split} \begin{array}{ll} T=\frac{2\,Z^{(r)} }{ Z^{(r)}+Z^{(l)}},& R=\frac{ Z^{(r)}-Z^{(l)} }{ Z^{(r)}+Z^{(l)}} \;. \\ \end{array} \end{split}\]

Let \(A_I\) denote the amplitude of the wave that hits the interface. The transmitted wave, with amplitude \(A_T = T A_I\), propagates to the right with wave speed \(c^{(r)}\). The reflected wave, with amplitude \(A_R = R A_I\), propagates to the left with wave speed \(c^{(l)}\).

Run Maxwell_1D_Interface.py. The numerical solution is integrated in time to \(t=1.8\) using the fourth order accurate Runge-Kutta method and time-step \(k=\frac{h}{10}\) to make the time-integration errors negligible. You can specify the number of grid points \(m\) and the order of accuracy of the SBP finite difference operator. You can also choose the value of \(\epsilon_r>0\) by changing the variable eps2 in the code. In the special case \(\epsilon_r=1\) we have the same medium on both sides of the interface. We recommend using at least \(m=101\) grid points.

Questions:

  1. What happens when \(\epsilon_r>1\)?

  2. What happens when \(\epsilon_r<1\) (but still positive)?

  3. What happens when you increase the number of grid intervals \(m\) for a given finite difference stencil? Does the numerical approximation improve?

  4. What happens when you increase the order of accuracy for a given \(m\)? Does the numerical approximation improve?

Coupling of two different media in 2D#

Consider Maxwell’s equations in 2D with the computational setup presented in the Figure below. The computational domain consists of two different media, where the computational domain is split into 4 adjacent blocks. This setup is solved using SBP finite differences. Let \(\mu_1=1\) and \(\epsilon_1=1\) denote the physical parameters in medium 1. In medium 2, the parameters are denoted \(\mu_2=1\) and \(\epsilon_2\). The computational domain is the square \([-2, 2] \times [-2, 2]\). The initial data are

\[ H^{(z)}=\exp{\left(-\left(\frac{x-x_0}{r_{*}}\right)^2-\left(\frac{y-y_0}{r_{*}}\right)^2\right)} \]

where \(r_*=0.1\), \(x_0=-0.5\) and \(y_0=0.5\). The electric field is initially set to zero.

../../_images/2D_Maxwell2.pdf

Here we will test 3 different sets of boundary conditions:

  1. Specifying the tangential electric field (i.e., either \(E^{(x)}\) or \(E^{(y)}\), depending on which boundary).

  2. Specifying the magnetic field \(H^{(z)}\).

  3. Specifying non-reflecting boundary conditions (particular linear combinations of electric and magnetic fields).

Start by running Maxwell_2D_Interface.py. The numerical solution is integrated in time to \(t=4\) using the fourth order accurate Runge-Kutta method and time-step \(k=\frac{h}{10}\). You can specify the number of grid points \(m\) and the order of accuracy of the SBP finite difference operator. You can also choose between the 3 types of boundary conditions listed above. You can also experiment with different values of \(\epsilon_2>0\) (eps2 in the code). A value \(\epsilon_2>1\) implies that the wave speed is lower in medium 2.

In the numerical tests we recommend using at least \(m=101\) (in each dimension and block, which corresponds to \(4\times m^2\) grid-points in total), but not more than roughly \(m=201\), since the runtime increases quickly with \(m\).

Try to answer the following questions:

  1. What happens to the solution when you choose non-reflecting boundary conditions? Compare this to the non-reflecting boundary conditions that you saw in 1D earlier.

  2. What happens when \(\epsilon_2>1\)?

  3. What happens when \(\epsilon_2<1\) (but still positive)?