# Interactive Educational Modules inScientific Computing

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(xc 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+1R ⁄ 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:

1. 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.
2. 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