Advection Equation
This module illustrates fully discrete finite difference methods for
numerically solving the advection equation. The advection equation
(also called the convection equation or one-way wave equation) in one
dimension is the partial differential equation
ut = −c
ux, where the solution
u(t, x) is a function of the time
variable t and the spatial variable x, and subscripts
indicate partial differentiation with respect to the given independent
variable. Considering the advection equation as a pure initial value
problem (a Cauchy problem) with initial time t = 0,
an appropriate initial condition is u(0, x) =
f(x) for a function f defined on the real
line. For this module the initial condition is taken to be a wave form
for which f(x) = 0 for x outside the
interval [−1, 1], and f is given by some
formula in the interval [−1, 1]. By the chain rule, an
exact solution to the advection equation is given by
u(t, x) = f(x − c
t); as time increases the initial condition retains its
shape while shifting to the right (or left, depending on the sign
of c) with velocity c. Computed numerical
solutions can be compared to this known exact solution to determine
their accuracy.
The user begins by selecting the constant velocity c and initial
function f. The initial function is plotted on the
interval [−3, 3]. Next, the user chooses the type
of finite difference method to be used. If the chosen type of method
uses a non-centered finite difference approximation for a spatial
derivative, the user can choose whether to use forward or backward
differences (either of which may be “upwind” or
“downwind”, depending on the sign of c). The user
may also choose between explicit and implicit versions of the method
type, if available. The stencil of the currently selected method is
shown below the method options, with the point being computed colored
red, the points it depends on colored blue, and the remaining points
colored black. Finally, the user specifies the step sizes in space and
time for the discrete mesh of points used in the finite difference
method. The stability limit shown below gives, for the chosen c and
space step size, the restriction on the size of the time step in order
for the numerical method to be stable. If the time step chosen
violates this restriction, then the numerical solution may oscillate
wildly and bear little resemblance to the true solution. Values of
infinity or zero for the stability limit indicate unconditional
stability or instability, respectively.
To view the numerical solution, the user chooses between
two-dimensional and three-dimensional display modes and then clicks
Start. The approximate solution is advanced time step by time
step, and the plot of the solution is updated accordingly. In
two-dimensional display mode, the solution at the current time is
plotted as a curve on the spatial interval [−3, 3],
and solution values at each new time step replace those at the previous
time step. The approximate solution is shown in green, and the exact
solution is shown in red for comparison. In three-dimensional display
mode, the approximate solution is plotted as a surface over the
space-time plane, and solution values at each new time step are added
to the existing graph, extending it forward along the time axis. The
solution continues to advance until time t = 1 is
reached or the user clicks Stop. When the solution process is
stopped before t = 1, it can be resumed by again clicking
Start. Clicking Reset clears any solution that may be
partially calculated and redisplays the initial condition, allowing the
user to select a different velocity constant, initial condition,
solution method, or step sizes.
Details:
All of the finite difference methods implemented here except the
one-sided trapezoid method are discussed in reference [2] below. The
one-sided trapezoid method uses the mean of the one-sided explicit
method and the one-sided implicit method. Thus, the scheme using
backward differences is (1 − R ⁄ 2)
ukn + R
⁄ 2 uk−1n
= (1 + R ⁄ 2)
ukn+1 −
R ⁄ 2
uk−1n+1, where
R = c Δt ⁄
Δx. The numerical methods are computed on the
interval [−5, 5], with only the interval
[−3, 3] displayed. For explicit methods, any
necessary boundary values are taken to be zero, based on the assumption
of a Cauchy problem and the limitations on the velocity constant. Many
of the implicit methods are implemented similarly. The one-sided
implicit methods are implemented with periodic boundary conditions to
avoid the possibility of a singular matrix when they are applied to
problems going the “wrong” direction.
References:
- Michael T. Heath, Scientific Computing,
An Introductory Survey, 2nd edition, McGraw-Hill, New York,
2002. See Section 11.1, especially Example 11.1 and Figure 11.1.
- J. W. Thomas, Numerical Partial Differential Equations: Finite
Difference Methods, Springer, New York, 1995. See Sections 5.3 and
5.4.
Developers: Evan VanderZee and Michael Heath