Numerical Analysis of an Advection Equation
This is a partial differential equation solved numerically using Matlab. There is a little bit of stability analysis and other stuff here, but my primary interest is in verifying all of the theory with a numerical scheme. I didn't go into a lot of detail about the analysis - but I will directly verify the results.
Advection is a process where a quantity (energy, mass, vorticity) moves from one region to another via the motion of a fluid like water or the atmosphere. Check out the plots below if advection is unclear: It is a process that is familiar to us in the regular world. Let's consider an advection equation with only one spatial dimension \( x\), velocity \( v\), and time \( t \) with some constant \( \alpha \), given by:
and initial boundary conditions given by
There is nothing special about the shape of the initial condition. I just picked it because it is easy to see it move on the interval \( [0, 10] \). There is no particular reason, either, to pick the interval \( [0, 10] \). The cool-looking function just fits well on that interval. Note that there is a periodic boundary condition. In other words the advection equation takes on the same values at \(x = 0 \) and \( x = 10 \).
We will run an Explicit Lax-Wendroff scheme.
the Explicit Lax-Wendroff scheme for \( u_t + \alpha u_x = 0 \) is
$$
\begin{align*}
u_j^{n+1} &= \frac{1}{2}\nu(1+\nu)u_{j-1}^n + (1 - \nu^2)u_j^n - \frac{1}{2}\nu(1 - \nu)u_{j+1}^n \\
&= u_j^n - \frac{\nu}{2}(u_{j+1}^n - u_{j-1}^n) + \frac{1}{2}\nu^2(u_{j+1}^n - 2u_j^n + u_{j-1}^n)
\end{align*}
$$
I made the substitution \( \displaystyle\nu = \alpha \frac{dt}{dx} \) for simplicity.
The scheme produces the stencil:
That scheme equation was derived directly using the 'forward-difference' or 'backward-difference' definition for partial differentiation. (That's the limit definition that freaks everyone out in calculus I.) The idea is to get an explicit expression for the next time step \( u_j^{n+1} \) from known quantities like \( u_j^n \) coming from the initial condition and step sizes in time and the spacial dimension. Doing this part is essential if we want to have any chance of getting a good solution for the problem. If we can't do it in terms of knowns we are possibly stuck with a bunch of unknowns - and just make more and more and more as we try to solve the problem!
This is a conditionally stable scheme - so it would be useful to know which values of \( t \), \( x \), and \( \alpha \) will produce numerically stable results. Here's the stability condition:
\begin{align*} u_j^{n+1} &= u_j^n - \frac{\nu}{2}\left( u_{j+1}^n - u_{j-1}^n \right) + \frac{\nu^2}{2}\left( u_{j+1}^n - 2u_j^n + u_{j-1}^n \right) \\ \lambda u_j^n &= u_j^n - \frac{\nu}{2} \left( u_j^n e^{ik(\bigtriangleup x)} - u_j^n e^{-ik(\bigtriangleup x)} \right) + \frac{\nu^2}{2}\left( u_j^n e^{ik(\bigtriangleup x)} - 2 u_j^n + u_j^n e^{-ik(\bigtriangleup x)} \right) \\ \lambda &= 1 - \frac{\nu}{2}\left( 2 i \sin(k \bigtriangleup x) \right) + \frac{\nu^2}{2}\left( 2 \cos(k \bigtriangleup x) - 2 \right) \\ \lambda &= 1 - \nu i \sin(k \bigtriangleup x) + \nu^2 (\cos(k \bigtriangleup x) - 1 ) \end{align*} We require \begin{align} \vert \alpha \vert \frac{\Delta \ t}{\Delta \ x} \leq 1 \end{align} This is important because the scheme is numerically stable if and only if the condition is satisfied. For example, if we pick too big a step size in time we will never be able to get a stable solution. It will blow up, as it were. Solutions will fly off to infinity... Dogs and cats will start sleeping together... A math teacher somewhere cries in class for no reason...We can write the scheme in matrix form. This is a very helpful thing to do here because we are going to use a computer to do all of the calculations... And the computer will be performing repeated multiplications using exactly this matrix. \begin{align*} \begin{bmatrix} u_1^{n+1} \\ \vdots \\ u_{j-1}^{n+1} \\ u_{j}^{n+1} \\ u_{j+1}^{n+1} \\ \vdots \\ u_{k}^{n+1} \end{bmatrix} = \begin{bmatrix} 1 - \nu^2 & \frac{1}{2}(\nu^2 - \nu) & & & & 0 & \frac{1}{2}(\nu^2 + \nu) \\ \frac{1}{2}(\nu^2 + \nu) & 1 - \nu^2 & \frac{1}{2}(\nu^2 - \nu) & & & 0 & 0 \\ 0 & \frac{1}{2}(\nu^2 + \nu) & \ddots & & & & \\ & & \frac{1}{2}(\nu^2 + \nu) & 1 - \nu^2 & \frac{1}{2}(\nu^2 - \nu) & & \\ & & & & \ddots & & \\ 0 & 0 & & & \frac{1}{2}(\nu^2 + \nu) & 1 - \nu^2 & \frac{1}{2}(\nu^2 - \nu) \\ \frac{1}{2}(\nu^2 - \nu) & 0 & & & & \frac{1}{2}(\nu^2 + \nu) & 1 - \nu^2 \end{bmatrix} \begin{bmatrix} u_1^{n} \\ \vdots \\ u_{j-1}^{n} \\ u_{j}^{n} \\ u_{j+1}^{n} \\ \vdots \\ u_{k}^{n} \end{bmatrix} \end{align*} The Matlab code more or less takes this matrix from math-land into Matlab-land and then pounds out solutions and plots everything.
Matlab Code: this is the main function for the Explicit Lax-Wendroff. The initial condition is a separate function so it can be reused.
function AdvectionStable(a,dt,dx,max_time)
% Creates function input x = (0 dx 2*dx ... 10)^T
% dx = .05;
x = 0:dx:10;
Nu = a*dt/dx;
n_dim = length(x);
time_array = 0:dt:max_time;
time_index = length(time_array);
% Find IC
[IC] = AdvectionStableIC(dx);
plot(x,IC)
pause(.1)
Ones = linspace(1,1,n_dim-1);
A = eye(n_dim) - Nu^2*eye(n_dim) + (.5)*(Nu^2 + Nu)*diag(Ones,-1) + (.5)*(Nu^2 - Nu)*diag(Ones,-1);
A(1,n_dim)=.5*(Nu^2+Nu);
A(n_dim,1)=.5*(Nu^2-Nu);
% Solved and plotted with Euler
U = IC';
for i = 0:1:time_index;
U = A*U;
plot(x,U(1:n_dim));
axis([0 10 -.1 2]);
pause(.01)
% Create snapshots to make a movie. The i+1 is necessary since we want to start M with the index of 1.
M(i+1) = getFrame(gcf);
end
% Produce 'AdvectionMovie.avi' using 'M' ** Remove comment marker to produce a movie **
% movie2avi(M,'AdvectionMovie.avi');
% For debugging
% dx
% dx
% a
% Nu
% A
end
function[IC] = AdvectionStableIC(dx)
% Creates function input x = (0 dx 2*dx ... 10)^T
x = 0:dx:10;
% length(x)
% IC = v(x,0)
IC = exp(-7*(x-2).^2) + exp(-(x-4).^2);
% To plot IC by itself:
plot(x,IC)
end
The code is pretty straightforward - and is fairly readable since I named the variables logical things like Nu \(= \nu \), etc. A few noteworthy things about Matlab: The built-in function 'eye(x)' creates an identity matrix. Also 'linspace(x,y,z)' can be used to create an evenly spaced vector. Both of these are useful shortcuts for constructing matrices and vectors - especially when you have to produce one of arbitrary size like this. The matrix and initial condition vector are created as soon as Matlab is told how big to make it. I always produce general code like this whenever possible. That also means that the function has to be called with all the necessary arguments.
function AdvectionStable(a,dt,dx,max_time)
This scheme produces very stable results depending upon the arguments.
Here is a movie produced using:AdvectionStable(1,0.05,0.05,10) AdvectionStable(0.5,0.05,0.05,10) AdvectionStable(2,0.05,0.05,10) AdvectionStable(1,0.055,0.05,10) I have written solutions for this PDE using different methods. I may add a page about that later.
Links to Matlab Code. Note that you must have both files for Matlab to run this code.
Main Lax-Wendroff Matlab Code
Initial Condition Lax-Wendroff Matlab Code
Depending on your browser you may have to 'save as.'
Some places to find out more:
- An excellent example of a movie produced by Matlab 'screenshots.' I used exactly this technique in the Matlab code.
- "Numerical Methods for the Solution of Partial Differential Equations" by Luciano Rezzolla (PDF)