Boundary conditions for the advection equation discretized by a finite difference method

Question

I am trying to find some resources to help explain how to choose boundary conditions when using finite difference methods to solve PDEs.

The books and notes which I currently have access to all say similar things:

The general rules governing stability in the presence of boundaries are far too complicated for an introductory text; they require sophisticated mathematical machinery

(A. Iserles A First Course in the Numerical Analysis of Differential Equations)

For example, when trying to implement the 2-step leapfrog method for the advection equation:

$u_i^{n+1} = u_i^{n-1} + \mu (u_{i+1}^n - u_{i-1}^n)$

using MATLAB

M = 100; N = 100;

mu = 0.5;

c = [mu 0 -mu];
f = @(x)(exp(-100*(x-0.5).^2));

u  = zeros (M, N);
x = 1/(M+1) * (1:M);

u(:,1) = f(x);
u(:,2) = f(x + mu/(M+1));

for i = 3:N
    hold off;
    u(:,i) = conv(u(:,i-1),c,'same') + u(:,i-2);
    plot(x, u(:,i));
    axis( [ 0 1 0 2] )
    drawnow;
end

The solution behaves nicely until it reaches the boundary, when it very suddenly starts behaving badly.

Where can I learn how to handle boundary conditions like this?

Answer 1

Sloede's response is very thorough and correct. I just wanted to add a few points to make it easier to grasp.

Basically, any wave equation has an inherent wave speed and direction. For a one-dimensional wave equation: $$ u_t + a u_x = 0 $$ the wave speed is the constant $a$ which determines not only the speed at which the information is propagating in the domain but also its direction. If $a>0$, information is going from left to right and if $a<0$ it's the other way around.

For the leap frog method, when you discretize the equations you get: $$\frac{u_i^{n} - u_i^{n-2}}{2 \Delta t} + a \frac{u_{i+1}^{n-1} - u_{i-1}^{n-1}}{2 \Delta x} = 0$$ or: $$u_i^{n} = u_i^{n-2} + \mu (u_{i+1}^{n-1} - u_{i-1}^{n-1}) $$ where $\mu = -a \Delta t/\Delta x$. In your case, $\mu > 0$ which translates to a wave going to the left. Now if you think about it, a wave that is traveling to the left, will only need a boundary condition at the right boundary since all the values to the left are being updated via their right neighbors. In fact specifying any value at the left boundary is inconsistent with the nature of the problem. In certain methods, like simple upwind, this is taken care of automatically since the scheme also involves only the correct neighbors in it's stencil. In other methods, like the leap frog, you have to specify some "correct" value.

This is usually done via extrapolation from the inside domain to find the missing value. For multi-dimensional and non-canonical problems, this involves finding all the eigen-vectors of the flux Jacobian to determine which parts of the boundary actually need boundary conditions and which parts require extrapolation.

Answer 2

General answer
Your problem is that you do not set (or even specify) the boundary conditions at all - your numerical problem is ill-defined.

Generally, there are two possible ways to specify the boundary conditions:

1. Set the boundary conditions by specifying $u_{0}$ and $u_{101}$ externally, e.g. through the exact solution.
2. Change the numerical stencil so it will use only interior information at the boundary.

Which way you go heavily depends on the physics of your problem. For wave-equation type problems one usually determines the eigenvalues of the flux Jacobian in order to decide whether external boundary conditions are needed, or whether the interior solution is to be used (this method is commonly called 'upwinding').

---

Quick fix
In your scheme, the numerical stencil needs the values of $u_{i-1}^n$ and $u_{i+1}^n$ to calculate a value at location $i$ for time $n+1$. Thus at the left domain boundary (i.e. $i = 1$) it would need the value of $u_{0}^n$, which you have not set in your code example. The same problem arises for $u_{100}^{n+1}$, which needs $u_{101}^n$ to build the stencil.

An easy way to get around would be to set $u_{1}^n$ and $u_{100}^n$ manually at each time step, by interpolating linearly from the interior solution. This is not really physical (as it does not take into account the actual wave speed), but will in your case lead to a first satisfactory solution behavior.

You can find a modified version of your source code below:

M = 100; N = 100;

mu = 0.5;

c = [mu 0 -mu];
f = @(x)(exp(-100*(x-0.5).^2));

u  = zeros (M, N);
x = 1/(M+1) * (1:M);

u(:,1) = f(x);
u(:,2) = f(x + mu/(M+1));

for i = 3:N
    hold off;
    %u(:,i) = conv(u(:,i-1),c,'same') + u(:,i-2);

    % Apply the numerical stencil to all interior points
    for j = 2:M-1
        u(j,i) = u(j,i-2) + mu*(u(j+1,i-1) - u(j-1,i-1));
    end

    % Set the boundary values by interpolating linearly from the interior
    u(1,i) = 2*u(2,i) - u(3,i);
    u(M,i) = 2*u(M-1,i) - u(M-2,i);

    plot(x, u(:,i));
    axis( [ 0 1 0 2] )
    drawnow;
end

Answer 3

So I've looked at this in some more detail, and it seems that this (at least in the basic cases which I'm handling) depends on the group velocity of the method.

The leapfrog method (for instance) is:

$u_{i}^{n+1} = u_{i}^{n-1} + \mu (u_{i+1}^n - u_{i-1}^n)$

Trying solutions of the form $u_k^n = e^{i( \zeta k \Delta x + \omega(\zeta) n \Delta t)}$ we find:

$e^{2 i \omega \Delta t} = 1 + \mu e^{i\omega \Delta t} ( e^{i \zeta \Delta x} - e^{-i \zeta \Delta x} )$

$\sin(\omega \Delta t) = \mu \sin (\zeta \Delta x)$

$\frac{d\omega}{d\zeta} = \frac{\cos (\zeta \Delta x)}{\sqrt{1-\mu^2 sin^2(\zeta \Delta x)}} \in [-1,1]$

Now we need to find out the group velocity of the boundary conditions:

My method: $u_1^{n+2} = u_1^n + \mu u_2^{n+1}$

We can compute the boundary group velocity as follows:

$2 i \sin(\omega \Delta t) = \mu e^{i \zeta \Delta x}$

so to find some group velocities which the boundaries allow we need to find:

$\omega = c \zeta$

$\cos(\zeta \Delta x) = 0, \mu \sin(\zeta \Delta x) = 2 \sin (\zeta c \Delta t)$

$\zeta = \frac{\pi}{2 \Delta x}$ would give $\mu = 2 \sin (\frac{c \mu \pi}{2})$ for which a solution for $c \in [-1,1]$ will exist. (For most choices of $\mu$ at least)

---

The solution which I've found in the literature is to take $u_0^{n+1} = u_1^{n}$ since this has a boundary wave number which lies outside $[-1,1]$.

I've still quite quite a bit more to read up about this before I understand it completely. I think the key words I'm looking for are GKS theory.

Source for all this A Iserles Part III notes

---

A clearer calculation of what I've done can be found here: http://people.maths.ox.ac.uk/trefethen/publication/PDF/1983_7.pdf

Answer 4

Guys I am very new to this site. Maybe this isn't the place to ask, but please forgive me as I am very new here :) I am having an extremely similar problem, the only difference being the starting function which, in my case, is a cosine wave. My code is this: clear all; clc; close all;

M = 1000; N = 2100;

mu = 0.5;

c = [mu 0 -mu]; f = @(x)1- cos(20*pi*x-0.025).^2; u = zeros (M, N); x=0:(1/M):0.05; u(1:length(x),1) = f(x); u(1:length(x),2) = f(x - mu/(M)); x=linspace(0,1,M);

for i = 3:N hold off;

% Apply the numerical stencil to all interior points
for j = 2:M-1
    u(j,i) = u(j,i-2) - mu*(u(j+1,i-1) - u(j-1,i-1));
end

% Set the boundary values by interpolating linearly from the interior
u(M,i) =  2*u(M-1,i-1) - u(M-2,i-1);

plot(x, u(:,i)); axis( [ 0 1.5 -0.5 2] ) drawnow; %pause end

There is already this code here, but for some reason, probably related to the cosine wave, my code fails :/ any help would be appreciated :) thanks!