Interactive Educational Modules in
Scientific Computing

Convection-Diffusion Equation

This module illustrates fully discrete finite difference methods for numerically solving the convection-diffusion equation. The convection-diffusion equation (sometimes called the transport equation) in one dimension is the partial differential equation ut + c ux = d uxx, where the solution u(t, x) is a function of the time variable t and the spatial variable x, subscripts indicate partial differentiation with respect to the given independent variable, and c and d are constants that control the balance between convection and diffusion. Considering this equation as an initial-boundary value problem with initial time t = 0, the solution u is defined for all x in a given interval [a, b] and any nonnegative t. An initial condition is given by a function f(x) defined on the interval [a, b], and the left and right boundary conditions are given by functions g(t) and h(t), respectively, each defined for any nonnegative t.

The user begins by selecting the convection constant c, the diffusion constant d, and an initial function. The initial function, represented by a truncated Fourier sine series approximation, is plotted on a predetermined interval. Depending on the selected constants, boundary conditions g(t) and h(t) are determined to define a problem with a known exact solution. Next the user chooses a fully discrete finite difference method to apply to the convection-diffusion equation. Methods not named after specific people use a two-word naming convention in which the first and second words describe the finite difference approximations to the time derivative ut and the convective term ux, respectively. The diffusive term uxx is approximated by centered differences in all cases. The stencil of the selected method is shown below, with the point being computed colored red, the other points used in the difference scheme colored blue, and any remaining points colored black. The red and blue points are drawn with three different levels of intensity to indicate how many of the terms in the convection-diffusion equation they are used to approximate. A point that is used in approximating only one term (e.g., the diffusive term uxx) is drawn with the lowest intensity, while a point used in approximating all three terms has the highest intensity. 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 constants c and d and the selected space step size, the restriction on the size of the time step in order for the specified 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. A value of infinity for the stability limit indicates unconditional stability.

To view the numerical solution, the user chooses between two-dimensional and three-dimensional display modes and 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 [a, b], 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 the time t = 1/4 is reached or the user clicks Stop. When the solution process is stopped before t = 1/4, 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 different parameters.

Details: All of the finite difference methods implemented are discussed in [1], in which each of the six is shown to be a special case of a general finite difference equation, depending on the choice of two weights, φ and θ. Unfortunately, that general equation (4.6) in [1] is incorrect. Taking the terms in the same order, the correct coefficients are θ (R + ψ), 2(1 + θ ψ), θ (Rψ), (1 − θ) (R + ψ), 2 (1 − ψ (1 − θ)), and −(1 − θ) (Rψ), where R = c Δt ⁄Δx and ψ = φ | R | + 2 d Δt ⁄ (Δx)2. Noye's optimal two-weight method should be implemented using the corrected equation (4.6), with weights given in equation (4.5.2), rather than from the specific finite difference equation equation (4.5.5) given in [1]. The first five methods in the menu are referred to in [1] as FTCS (forward-time, centered-space), US (upstream), LW (Lax-Wendroff), BTCS (backward-time, centered-space), and CNT (Crank-Nicolson type), respectively.

References

  1. John Noye, Finite difference methods for solving the one-dimensional transport equation, in Numerical Modeling: Applications to Marine Systems, edited by J. Noye, Elsevier, 1987, pages 231-256.
  2. John C. Strikwerda, Finite Difference Schemes and Partial Differential Equations, 2nd edition, SIAM, Philadelphia, 2004. See Section 6.4.

Developers: Evan VanderZee and Michael Heath