Homework

Homework #

PDE Solvers: Finite Difference Methods 2 #

Table of Contents #

%matplotlib inline
import numpy as np
import scipy.linalg as sl
# we'll start using sparse matrices here
import scipy.sparse as sp
# and linear algebra functions designed for sparse matrices
import scipy.sparse.linalg as spla
import matplotlib.pyplot as plt
# the following allows us to plot triangles indicating convergence order
from mpltools import annotation

# font sizes for plots
plt.rcParams['font.size'] = 12
plt.rcParams['font.family'] = 'sans-serif'
plt.rcParams['font.sans-serif'] = ['Arial', 'Dejavu Sans']

Homework - Upwind vs central discretrisation #

Use our codes BVP_AD_central and BVP_AD_upwind to compute and plot error convergence for these two schemes for our test problem from the lecture with a Peclet number of 10.

Hence, establish that even though the upwind scheme is coupled with a second-order scheme for diffusion, it is the first-order contribution to the overall error which ultimately dominates and so the overall scheme is only first-order accurate.

Homework - Larger stencils / higher orders of accuracy [\(\star\)] [from Comp Math]#

We saw in the lecture first- and second-order approximations to the first derivative.

By making use of more points (expanding the stencil) it is possible to come up with approximations of arbitrary order.

Tables of coefficients, assuming a uniform mesh, can be found at: https://en.wikipedia.org/wiki/Finite_difference_coefficient

Try expanding our example from the lecture where we computed the first derivative of \(\sin(x)\) at the location \(x=0.8\) and plotted the error as a function of \(\Delta x\), with some other examples of finite difference stencils from the tables at that web page.

[NB. the part of the question above is quite simple - in the sample solution I demonstrate an implementation of the third order accurate scheme; the part of the question below is more complicated].

Note that in the section of that web page titled “Arbitrary stencil points” it also gives a matrix system which can be solved for the finite difference coefficients which provides an approximation of arbitrary order derivatices on an arbitrary stencil.

For a more in-depth description and derivation take a look at: http://web.media.mit.edu/~crtaylor/calculator.html (click on “How does it work?”).

Write some code to construct and solve this matrix system, and hence extend the convergence plot we presented above for the first derivative of \(\sin(x)\) and \(x=0.8\) for orders 1-6.

Homework - Stability study for advection-diffusion using FTCS #

Play around with the example from the lecture in the section “Stability study for advection-diffusion”. In particular play around with the physical and numerical parameters and investigate further the interplay between these, the locations of eigenvalues in the complex plane, and the observed stability of the numerical solver.

Homework - An analytical solution to advection-diffusion [from Comp Math]#

An exact solution to the advection-diffusion equation is given by

\[ c(x,t) = \frac{1}{\sqrt{4\pi \, \kappa \,t}}\exp\left (-\frac{(x-Ut)^2}{4\kappa \, t}\right).\]

Note that this holds in an infinite domain - we can make use of it only as long as our numerical solution is far away from boundaries, or in a periodic domain as long as the solution behaviour does not start to encroach on itself.
The initial condition that this solution corresponds to is a Dirac-delta function. We clearly can’t represent this on our mesh as it has infinitesimal thickness - all sorts of problems would arise if we tried to. But what we can do is assume that our simulation starts at some \(t>0\) and initialise our problem with the corresponding exact solution (which for \(t>0\) is a Gaussian and for sufficiently large times can be represented on our mesh).

Write a function to evaluate this exact analytical solution and use it to initialise a simulation and plot a comparison between the exact and analytical solution at later times.

I suggest you use periodic boundary conditions (e.g. your starting point could be our function from the lecture solve_adv_diff_Gaussian2).

First try selecting some appropriate physical and numerical parameters such that you have a stable solution and avoid boundaries. Then try refining the spatial and temporal mesh to see if your solution gets closer to the exact solution.

Homework - Pure advection - impact of time step [from Comp Math]#

Write a solver for a pure advection problem in a periodic domain using forward Euler with upwind differences in space and with a Gaussian initial condition.

First test that it works/fails based on the CFL condition.

Then update the solver such that it selects the timestep automatically to satisfy a user-defined target Courant number, e.g. look to enforce a Courant number of 0.8.