Numerical Boundary Conditions #
The linear advection equation is
\[\pdv{u}{t} + a \pdv{u}{x} = 0\]For it to be a well-posed mathematical problem, you also need the Cauchy data (boundary and initial conditions). From the order of the differential operators, we know that we need one boundary condition (constraint on \( u(x_b) \) for all \( t \) ) and one initial condition (constraint on \( u(x) \) at \( t = 0 \)).
If we are solving the heat diffusion equation
\[\pdv{T}{t} = \alpha \pdv{^2 T}{x^2}\]Then we need 1 initial condition and 2 boundary conditions.
Back to the linear advection equation, if we say
\[u(t = 0, x) = f(x) \quad \text{and} \quad u(t, x=0) = g_D(t) \qquad \text{(Dirichlet BC)}\]For wave speeds less than 0 (\( a < 0 \)), then we would instead specify \( u(t, x=L) \).
The way boundary conditions are implemented in the Forward Euler algorithm (as an example):
\[u_j ^{n+1} = u_j ^n - \frac{ a \Delta t}{\Delta x} (u_j ^n - u_{j-1} ^n)\]- We begin by specifying \( u_j ^0 = f(x_j) \), i.e. we populate the solution at \( t=0 \).
- Sweep through the domain starting at \( j=1 \) towards \( j=J \)
where
\[u_0 ^n = g_D(t^n)\]- Notice, at \( j=J \)
There is no boundary condition required since there is no data required from outside of our domain. If we do require such data, we’ll need to do something special here.
Leap-Frog BC #
Recall our Leap-frog algorithm
\[u_j ^{n+1} = u_j ^{n-1} - \frac{a \Delta t}{\Delta x} (u_{j+1} ^n - u_{j-1} ^n)\]With assumed given initial condition
\[u_j ^0 = f(x) = u(x, t=0)\]What happens on the very first step from \( n=0 \) to \( n=1 \)?
\[u_j ^1 = u_j ^{-1} - \frac{a \Delta t}{\Delta x} (u_{j+1} ^0 - u_{j-1} ^0)\]We’ve got a starting problem: we need to know \( u_j ^{-1} \), which is outside of our problem domain. In order to take the first time step, we bootstrap the algorithm by any some other algorithm which is not centered in time. Generally we use Lax-Wendroff to give us an additional initial condition \( u_j ^1 \). With both \( u_j ^0 \) and \( u_j ^1 \), we can proceed with any second-order time-centered algorithm.
How about the boundary conditions near the edges of the domain?
\[u_1 ^{n+1} = u_1 ^{n-1} - \frac{a \Delta t}{\Delta x} (u_2 ^n - u_0 ^n) \qquad (x = 0)\] \[u_J ^{n+1} = u_J ^{n-1} - \frac{a \Delta t}{\Delta x}(u_{J+1} ^n - u_{J-1} ^n) \qquad (x = L)\]We’re now being asked for information outside of the problem domain \( u_0 ^n \) and \( u_{J+1} ^n \). This goes back to the modified PDE that we’re actually solving. The linear advection equation is first-order in \( x \) so we only need a single boundary condition to solve it, but the modified PDE that we are solving with Leap-frog has a second-order diffusive term, so we require an additional boundary condition.
Our solution to the boundary problem is called a Numerical Boundary Scheme. As a zeroth-order scheme, we can just extend our solution at the boundary out one grid point
\[(j = J)\qquad u_{J+1} ^n = u_J \rightarrow u_J ^{n+1} - \frac{ a \Delta t}{\Delta x} (u_J ^n - u_{J-1} ^n)\]We could also make use of a first-order extrapolation
\[u_{J+1} ^n = 2 u_J ^n - u_{J-1} ^n \rightarrow u_J ^{n+1} = u_J ^{n-1} - \frac{2 a \Delta t}{\Delta x} (u_J ^n - u_{J-1} ^n)\]Another method that’s particularly well suited to leap-frog is
\[u_{J+1} ^n = u_J ^{n-1}\]The advantage of this method is that in a centered scheme it preserves information flow. It assumes that information from \( (n-1, J) \) will flow to \( (n, J+1) \), and we can use that information to compute \( (n+1, J) \). This helps to prevent unphysical reflections at the boundary.
Another concept which addresses unphysical boundary effects that we will cover later is called PML (Perfectly Matched Layer).
NBS are needed in all algorithms which use central difference operators.
Neumann Boundary Conditions #
Instead of specifying the value of \( u \) at the boundary, we can specify the value of \( \pdv{u}{x} \) instead
\[\left.\pdv{u}{x}\right|_{x=0} = g_N(t)\] \[\rightarrow \frac{u_1 ^n - u_0 ^n}{\Delta x} = g_N(t) \rightarrow u_0 ^n = u_1 ^n - \Delta x g_N(t)\]What happens when we put this into the Forward Euler algorithm
\[u_1 ^{n+1} = u_1 ^n - \frac{a \Delta t}{\Delta x} \left(u_1 ^n - \left[u_1 ^n - \Delta x g_N(t)\right]\right) = u_1 ^n - a \Delta t g_N (t)\]So the spatial difference has been entirely replaced by the boundary condition, which is about what we would expect. It means that we can essentially remove \( j=1 \) from our algorithm, since we can write down the value of \( u_1 ^n \) for all \( n \) only using the BC and IC.
Numerical BC in Implicit Algorithms #
The way we implement boundary conditions in implicit algorithms is by incorporating the BC into the operator matrix. Looking at the \( \theta \)-method:
\[- \theta \frac{a \Delta t}{\Delta x} u_{j-1} ^{n+1} + u_j ^{n+1} + \theta \frac{a \Delta t}{\Delta x} u_{j+1} ^{n+1} = u_j ^n (1 - \theta) \frac{a \Delta t}{2 \Delta x} (u_{j+1} ^n - u_{j-1} ^n)\]Consider \( j=1 \), and define \( \nu \equiv \frac{a \Delta t}{2 \Delta x} \) to clean things up
\[(j = 1) \quad u_1 ^{n+1} + \theta \nu u_2 ^{n+1} = u_1 ^n - (1 - \theta) \nu u_2 ^n + \nu [\theta u_0 ^{n+1} + (1 - \theta) u_0 ^n]\]When we have Dirichlet boundary conditions, the terms in square brackets \( u_0 ^{n+1}, u_0 ^n \) are the ones we would replace by our boundary condition. Writing out the operator matrix for the algorithm, the boundary condition is just an additional term we add
\[\begin{bmatrix} 1 & \theta & \nu & 0 & 0 & \ldots \\ -\theta \nu & 1 & \theta \nu & 0 & \ldots \\ \ldots \end{bmatrix} \begin{bmatrix} u_1 ^{n+1} \\ u_2 ^{n+1} \\ \ldots \end{bmatrix} \] \[= \begin{bmatrix} 1 & -(1-\theta) \nu & 0 & \ldots & \ldots \\ (1-\theta) \nu & 1 & -(1-\theta) \nu & 0 \ldots \\ \ldots \end{bmatrix} \begin{bmatrix} u_1 ^n \\ u_2 ^n \\ \ldots \end{bmatrix} + \begin{bmatrix} \nu \theta g_D(t^{n+1}) + \nu (1 - \theta ) g(t^n) \\ 0 \\ 0 \\ \ldots \end{bmatrix}\]Note that now \( u_0 \) does not appear in the solution vector, so whenever we need to go and plot the solution, we need to add the initial condition back in.
Now consider the Neumann BC for the \( \theta \)-algorithm.
\[u_0^n = u_1 ^n - \Delta x g_N (t)\]At \( j=1 \)
\[\theta \frac{a \Delta t}{\Delta x} (\Delta x g_N(t^{n+1})- \theta\frac{ a \Delta t}{2 \Delta x} u_1 ^{n+1}+ u_1 ^{n+1} + \theta \frac{ a \Delta t}{2 \Delta x} u_2 ^{n+1} = u_1 ^n - (1 - \theta) \frac{a \Delta t}{2 \Delta x} (u_2 ^n - u_1 ^n + \Delta x g_N (t^n))\]Notice that
\[\theta \frac{a \Delta t}{\Delta x} (\Delta x g_N(t^{n+1})- \theta\frac{ a \Delta t}{2 \Delta x} u_1 ^{n+1} = - \theta \frac{ a \Delta t}{2\Delta x} u_0 ^{n+1}\]So if we look at our operator matrix,
\[\begin{bmatrix} 1 - \theta \nu & \theta \nu & 0 & 0 \ldots \\ - \theta \nu & 1 & \theta \nu & 0 \ldots \\ \ldots \end{bmatrix} \begin{bmatrix} u_1 ^{n+1} \\ u_2 ^{n+1}\\ \ldots \end{bmatrix} \begin{bmatrix} u_1 ^{n+1} \\ u_2 ^{n+1} \\ \ldots \end{bmatrix}\] \[= \begin{bmatrix} 1 + (1 - \theta)\nu & -(1-\theta) \nu & 0 & \ldots \\ (1 - \theta) \nu & 1 & - (1 - \theta) \nu & 0 \ldots \\ \ldots \end{bmatrix} + \begin{bmatrix} + - \theta \nu \Delta x g_N(t^{n+1}) - (1 - \theta) \nu \Delta x g_N(t^n) \\ 0 \\ \ldots \end{bmatrix}\]In comparison with Dirichlet BC, now we’ve actually changed both the operator matrix and the inhomogeneity vector (instead of just adding a term to the inhomogeneity vector).
Since the \( \theta \)-method uses central difference operators, NBS is also required at the other end of the domain at \( j=J \). We have the same options to choose from (0-th order or 1st-order extrapolation).
\[u_{J+1} ^n = u_J ^n \qquad \text{0th Order}\] \[\rightarrow - \theta \nu u_{J-1} ^{n+1} + ( 1 + \theta \nu ) u_J^{n+1} = u_J ^n - (1 - \theta) \nu( u_J ^n - u_{J-1} ^n)\] \[A = \begin{bmatrix} \ldots & \ldots & \ldots & \ldots & \ldots \\ \ldots & 0 - \theta & \nu & 1 & \theta \nu \\ \ldots & \ldots & 0 & - \theta \nu & 1 + \theta nu \end{bmatrix}\]There is again no row for \( u_{J+1} \). How about 1st order:
\[u_{J+1} ^n = 2 u_J ^n - u_{J-1}^n\] \[\rightarrow -2 \theta \nu u_{J-1} ^{n+1} + (1 + 2 \theta \nu ) u_J ^{n+1} = u_J ^n - 2 (1 - \theta) \nu (u_J ^n - u_{J-1} ^n) \] \[A = \begin{bmatrix} \ldots & \ldots & \ldots & \ldots \\ \ldots & - \theta \nu & 1 & \theta \nu \\ \ldots & 0 & - 2 \theta \nu & 1 + 2 \theta \nu \end{bmatrix}\]Periodic Boundary Conditions #
Periodic boundary conditions alter the operator matrix in a fairly straightforward manner. We simply set
\[u_0 = u_J \qquad u_{J+1} = u_1\] \[A = \begin{bmatrix} 1 & \theta \nu & 0 & \ldots & \ldots & \ldots & 0 & - \theta \nu \\ -\theta \nu & 1 & \theta \nu & 0 & \ldots& \ldots & \ldots & 0 \\ \ldots & \ldots& \ldots& \ldots \\ 0 & \ldots & \ldots & \ldots & \ldots &- \theta \nu & 1 & \theta \nu \\ \theta \nu & 0 & \ldots & \ldots & \ldots & 0 & - \theta \nu & 1 \end{bmatrix}\]