Skip to content

Numerical Boundary Conditions

The linear advection equation is

∂u∂t+a∂u∂x=0 \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(xb) u(x_b) for all t t ) and one initial condition (constraint on u(x) u(x) at t=0 t = 0 ).

If we are solving the heat diffusion equation

∂T∂t=α∂2T∂x2 \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)andu(t,x=0)=gD(t)(Dirichlet BC) 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 a < 0 ), then we would instead specify u(t,x=L) u(t, x=L) .

The way boundary conditions are implemented in the Forward Euler algorithm (as an example):

ujn+1=ujn−aΔtΔx(ujn−uj−1n) u_j ^{n+1} = u_j ^n - \frac{ a \Delta t}{\Delta x} (u_j ^n - u_{j-1} ^n)

  1. We begin by specifying uj0=f(xj) u_j ^0 = f(x_j) , i.e. we populate the solution at t=0 t=0 .
  2. Sweep through the domain starting at j=1 j=1 towards j=J j=J

    u1n+1=u1n−aΔtΔx(u1n−u0n) u_1 ^{n+1} = u_1 ^n - \frac{a \Delta t}{\Delta x} (u_1 ^n - u_0 ^n)

    where

    u0n=gD(tn) u_0 ^n = g_D(t^n)

  3. Notice, at j=J j=J

    uJn+1=uJn−aΔtΔx(uJn−uJ−1n) u_J ^{n+1} = u_J ^n - \frac{a \Delta t}{\Delta x} (u_J ^n - u_{J-1} ^n)

    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

ujn+1=ujn−1−aΔtΔx(uj+1n−uj−1n) 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

uj0=f(x)=u(x,t=0) u_j ^0 = f(x) = u(x, t=0)

What happens on the very first step from n=0 n=0 to n=1 n=1 ?

uj1=uj−1−aΔtΔx(uj+10−uj−10) 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 uj−1 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 uj1 u_j ^1 . With both uj0 u_j ^0 and uj1 u_j ^1 , we can proceed with any second-order time-centered algorithm.

How about the boundary conditions near the edges of the domain?

u1n+1=u1n−1−aΔtΔx(u2n−u0n)(x=0) u_1 ^{n+1} = u_1 ^{n-1} - \frac{a \Delta t}{\Delta x} (u_2 ^n - u_0 ^n) \qquad (x = 0)

uJn+1=uJn−1−aΔtΔx(uJ+1n−uJ−1n)(x=L) 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 u0n u_0 ^n and uJ+1n 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 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)uJ+1n=uJ→uJn+1−aΔtΔx(uJn−uJ−1n) (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

uJ+1n=2uJn−uJ−1n→uJn+1=uJn−1−2aΔtΔx(uJn−uJ−1n) 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

uJ+1n=uJn−1 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) (n-1, J) will flow to (n,J+1) (n, J+1) , and we can use that information to compute (n+1,J) (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 u at the boundary, we can specify the value of ∂u∂x \pdv{u}{x} instead

∂u∂x∣x=0=gN(t) \left.\pdv{u}{x}\right|_{x=0} = g_N(t)

→u1n−u0nΔx=gN(t)→u0n=u1n−ΔxgN(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

u1n+1=u1n−aΔtΔx(u1n−[u1n−ΔxgN(t)])=u1n−aΔtgN(t) 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 j=1 from our algorithm, since we can write down the value of u1n u_1 ^n for all n 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:

−θaΔtΔxuj−1n+1+ujn+1+θaΔtΔxuj+1n+1=ujn(1−θ)aΔt2Δx(uj+1n−uj−1n) - \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 j=1 , and define ν≡aΔt2Δx \nu \equiv \frac{a \Delta t}{2 \Delta x} to clean things up

(j=1)u1n+1+θνu2n+1=u1n−(1−θ)νu2n+ν[θu0n+1+(1−θ)u0n] (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 u0n+1,u0n 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

[1θν00…−θν1θν0……][u1n+1u2n+1…] \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}

=[1−(1−θ)ν0……(1−θ)ν1−(1−θ)ν0……][u1nu2n…]+[νθgD(tn+1)+ν(1−θ)g(tn)00…] = \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 u0 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.

u0n=u1n−ΔxgN(t) u_0^n = u_1 ^n - \Delta x g_N (t)

At j=1 j=1

θaΔtΔx(ΔxgN(tn+1)−θaΔt2Δxu1n+1+u1n+1+θaΔt2Δxu2n+1=u1n−(1−θ)aΔt2Δx(u2n−u1n+ΔxgN(tn)) \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

θaΔtΔx(ΔxgN(tn+1)−θaΔt2Δxu1n+1=−θaΔt2Δxu0n+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} = - \theta \frac{ a \Delta t}{2\Delta x} u_0 ^{n+1}

So if we look at our operator matrix,

[1−θνθν00…−θν1θν0……][u1n+1u2n+1…][u1n+1u2n+1…] \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}

=[1+(1−θ)ν−(1−θ)ν0…(1−θ)ν1−(1−θ)ν0……]+[+−θνΔxgN(tn+1)−(1−θ)νΔxgN(tn)0…] = \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 j=J . We have the same options to choose from (0-th order or 1st-order extrapolation).

uJ+1n=uJn0th Order u_{J+1} ^n = u_J ^n \qquad \text{0th Order}

→−θνuJ−1n+1+(1+θν)uJn+1=uJn−(1−θ)ν(uJn−uJ−1n) \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=[………………0−θν1θν……0−θν1+θnu] 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 uJ+1 u_{J+1} . How about 1st order:

uJ+1n=2uJn−uJ−1n u_{J+1} ^n = 2 u_J ^n - u_{J-1}^n

→−2θνuJ−1n+1+(1+2θν)uJn+1=uJn−2(1−θ)ν(uJn−uJ−1n) \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=[……………−θν1θν…0−2θν1+2θν] 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

u0=uJuJ+1=u1 u_0 = u_J \qquad u_{J+1} = u_1

A=[1θν0………0−θν−θν1θν0………0…………0…………−θν1θνθν0………0−θν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}