Complex Boundary Value Problems of Nonlinear Differential Equations: Theory, Computational Methods, and ApplicationsView this Special Issue
General Formulation of Second-Order Semi-Lagrangian Methods for Convection-Diffusion Problems
The general formulation of the second-order semi-Lagrangian methods was presented for convection-dominated diffusion problems. In view of the method of lines, this formulation is in a sufficiently general fashion as to include two-step backward difference formula and Crank-Nicolson type semi-Lagrangian schemes as particular ones. And it is easy to be extended to higher-order schemes. We show that it maintains second-order accuracy even if the involved numerical characteristic lines are first-order accurate. The relationship between semi-Lagrangian methods and the modified method of characteristic is also addressed. Finally convergence properties of the semi-Lagrangian finite difference schemes are tested.
There have been many numerical methods developed to deal with convection-dominated diffusion problems, and among them characteristic-based methods are the most popular. Numerical methods that follow characteristic lines backwards in time and then interpolate at their feet date back to the work of Courant et al. . As one kind of characteristics-based method, the semi-Lagrangian (SL) method was introduced in the beginning of the 1980s by Robert . Its basic idea is to discretize the Lagrangian derivative of the solution instead of the Eulerian derivative. This technique can increase significantly the maximum allowable time step while maintaining the efficiency of symmetric solvers. The SL methods have been extensively applied in numerical simulations of models for weather forecast and oceanography (see Staniforth and Côté  and Smolarkiewicz and Pudykiewicz  for review). As another kind of characteristic-based methods, the modified method of characteristic (MMOC) was introduced by Douglas Jr. and Russell  at roughly the same time as the SL method and has been extensively implemented in numerical simulations of fluid flows in porous media (see [6–8]) and many other transport problems (see [9–12] for review).
Strang  pointed out that the first-order methods were often too crude and the third-order methods too complicated. The computations are thus made expensive either by the fine mesh required by a first order scheme in order to provide enough detail or else by the delicate differencing which maintains a high-order accuracy. Second-order schemes are the obvious compromise.
Based on MMOC, in the beginning of 1980s Ewing and Russell  introduced the backward difference formula (BDF) of characteristic for linear constant-coefficient convection-diffusion problems. Afterwards characteristic schemes of Crank-Nicolson (CN) type for convection-diffusion equations and the Navier-Stokes (NS) equations were studied by Rui and Tabata  and Notsu and Tabata , respectively. The SL-BDF schemes for incompressible NS problems were proposed by Boukir et al. [17, 18]. The SL-CN methods were proposed for convection-diffusion equations and/or the incompressible NS equations by Bermúdez et al. [19, 20], Fourestey and Piperno , Xiu and Karniadakis , and Xiu et al. . Al-Lawatia et al.  and Falcone and Ferretti  presented and analyzed the single-step high-order semi-Lagrangian schemes of the Runge-Kutta type. More research was given by Bermejo and Conde , Xiao and Yabe , and Toda et al. .
Up to now we have not seen the SL method to be formulated in a sufficiently general fashion. In this paper we present a general second-order SL formulation for convection-dominated diffusion problems. By the method of lines approach, this formulation includes most of the second-order SL schemes mentioned previously. Extension of it to higher-order schemes is also addressed. Our development of the general SL formulation is mainly motivated by the work of Rui and Tabata  and Notsu and Tabata . They treated convection-diffusion and the NS problem, respectively, by MMOC. They emphasized that an “additional correction term” was indispensable to maintain the second-order accuracy of the MMOC schemes. We will see that the correction term of the MMOC schemes is actually a natural term of the SL schemes.
First, we show that by the method of lines (MOL) approach, the general SL formulation includes the SL-BDF2 and the SL-CN schemes as specific ones. Second, we verify that the formulation maintains second-order accuracy even if the involved numerical characteristic lines are first-order accurate. Third, we show that MMOC can be considered as a special version of the SL method. Finally, combining finite difference discretization in spaces, a fully discretized SL scheme is presented. Numerical tests demonstrate that the SL finite difference scheme is second-order convergent.
The outline of the rest of this paper is as follows. In Section 2 the forming of the SL schemes is recalled and the relationship between the first-order SL methods and MMOC is addressed. In Section 3 the general SL formulation is derived and its second-order consistency is proved. The relation of the second-order SL method and MMOC is also addressed. In Section 4 the finite difference SL schemes are derived. In Section 5 several SL finite difference schemes are applied to the Gaussian hill problem. In Section 6 conclusions of this paper are drawn.
2. The SL Methods and the MMOC
In this section we recall some details of the construction of the SL method and MMOC. Also we study the relationship between them.
Let be a bounded domain in with Lipschitz boundary , and let be a positive constant. Without loss of generality, we consider the initial boundary value problem: find such thatwhere is a positive constant, , and are given functions. For simplicity, we assume that is a divergence-free velocity field, that is, , and vanishes on the boundary .
Characteristic line is the trajectory of a fluid particle. The travel of the particle is associated with the velocity field . For a given point , the characteristic line through is represented by the vector function , which solves the initial value problemHere is the position of the particle on the characteristic line at time . The particle is located at at time . We assume that is Lipschitz continuous on with respect to the first variable. By the theory of ODEs, the characteristic line is well defined. By the chain rule we have where is the composition of functions with respect to the first argument of . Similar to the deriving of and in , (1a) can be written to Lagrangian form
For a positive integer , let be time step length and , for . Let denote the characteristic line on (or for two-step methods) through . Thus (4) can locally be written to
Tracking the particle backward from to along the characteristic line on (or from to on for two-step methods) corresponds to the backward solution of the following Cauchy problem: For function , denote and for . By applying the backward Euler’s method to (5), it follows that where represents the approximate of . Note that in practice, (6) can usually be solved approximately. By applying the backward Euler’s method to (6), it follows that and (similarly we see that for the two-step methods). In (7) with being replaced by , we have
On the other hand, (9) can be derived by MMOC . In fact, with , let denote the direction vector , and define the operator with . So (1a) can be written to the form Using the backward difference quotient, we have Substituting (12) into (11), we obtain (9). Thus MMOC can be considered as a special version of SL. In the next section we will see that similar relation exists for second-order case.
3. General Second-Order SL Formulation
In this section we present a general second-order SL formulation and show that by MOL this formulation includes several widely used schemes as specific ones.
Let us first consider the abstract ODEs of the following form: given the Hilbert space and , find such that where . A general second-order scheme for (13) can be of the form (see [29, 30]) As usual we denote scheme (14) as , where and are the characteristic polynomials of the scheme, with We further assume that is normalized by and satisfies the following conditions: for , , such that is Dahlquist and -stable .
For a fixed and , let . Then the scheme for (5) is of the form With being replaced by the approximate characteristic line , we have the analogue of (17): Scheme (18) is the general second-order SL formulation which by MOL includes all the previously introduced second-order SL schemes. Next we will prove that it includes the SL-BDF2 schemes and the SL-CN schemes as specific cases.
3.1. The BDF2 Scheme
In (14), let ; then from (16), we have , , . Substituting these coefficients into (14), we obtain the BDF2 scheme: Due to its stability and damping properties, (19) is one of the most popular second-order schemes . Substituting the previously coefficients into (17), we get the SL analogue of (19): Furthermore, in (20) with and being, respectively, approximated by , , we have This is just the multistep characteristic scheme derived based on MMOC in . In fact, using MMOC (analogous to (12)) we see that Substituting (22) into (11), we obtain (21).
The SL-BDF2 schemes were presented and analyzed in [14, 17, 18, 21]. It is deserved to note that though the involved approximate characteristic line is first-order accurate, the resulting SL-BDF2 scheme (21) maintains second-order accurate.
3.2. The CN Scheme
In (14), let From (14) we obtain the CN scheme From (17) we get the SL analogue of (24): and with being replaced by , we have
Remark 1. It is deserved to note that Rui and Tabata  and Notsu and Tabata  called -term in (26) the “additional corrected term,” since it is introduced from “outside” to recover the second-order consistency of the MMOC schemes. But in view of the previous discussion, the “additional corrected term” is a natural one in the SL schemes. Thus we think the SL method is more general than MMOC.
3.3. The Consistency of the General Formulation
Now we show that both (17) and (18) with first-order approximate characteristic lines are second-order accurate. Let us denote for . By the theory of ODEs (see ), using the Taylor’s expansion, from (13) and (14) we have If the following conditions hold then . It is easy to see that conditions (29) and (16) are equivalent if . In (13), let From (28)–(30) it follows that where is the exact solution of (5). Thus we have confirmed the following proposition.
Proposition 2. Suppose that and in (1a)–(1c) are smooth functions in , and is the exact characteristic line which solves (6). Then (17) is second-order consistent with (5).
Analogous to Proposition 2, if first-order characteristic lines are involved, then the following proposition holds.
Proposition 3. Suppose that and in (1a)–(1c) are smooth functions in . In (18), if , , then the scheme is second-order consistent with (5).
Proof. Analogous to the deriving of and in , for (5) can be written to the form
With being replaced by and being the approximate of , we obtain -perturbation of (33) of the form
Rewrite (36) to the Lagrangian form
Thus (37) is the second-order approximation to (5).
Corresponding to (30), let Similar to the proof of Proposition 2, we can see that (32) is second-order consistent with (37). We deduce that (32) is second-order consistent with (5).
Remark 4. In (32) higher-order numerical characteristic lines are usually preferred, though the characteristic line computed by the backward Euler’s method can retain the second-order accuracy. For example, Bermúdez et al. [19, 20], Rui and Tabata , and Notsu and Tabata  computed the numerical characteristic by the higher-order Runge-Kutta methods. In the following numerical tests we will see that the first-order characteristic line is too coarse to ensure reasonable convergence of the SL-CN schemes.
4. The SL Finite Difference Method
In this section we present a full-discretized SL formulation which combines finite difference for spatial discretizations. We also numerically verify the convergence of the formulation.
First, we build the finite difference scheme for (5). Assume that is a unit square, with boundary . Denote , with For the partition of , we denote and . Let , with Let , with Let , with . By the transformation given in Appendix section, term in (5) changes to expression that consists of , , , , , , , , , , and so forth. If , are polynomials of of degrees not more than one, then it holds that (see Appendix) Using centered difference to discretize these partial derivatives, we obtain the finite difference approximation of (5) as follows: find such that where is the approximate of . Applying (17) to (43), with representing the approximate of , we obtain the second-order-in-time finite difference scheme: In (44) with being replaced by , representing the approximate of , we finally discretize (1a)–(1c) as follows: The computation of the starting values of this multistep scheme is similar to the general Eulerian schemes. Since and in (45) are not grid points of , interpolations are needed. In the following numerical tests, we will use cubic interpolations (CI) and cubic spline interpolations, respectively (CSI).
5. Numerical Results
In this section, we test some specific cases of (45). The problem of rotating Gaussian hill has been widely used to test numerical schemes for convection-diffusion equations. In (1a)–(1c), let , , , , and or . The Dirichlet boundary conditions and initial condition of (1a)–(1c) are given such that the exact solution is where , , and . The exact characteristic is given by By the Euler’s method the approximate characteristic is given by
Note that since the velocity does not satisfy the nonflow boundary condition (1b), the characteristics and its approximations are not necessarily contained in . We assume that outside of (just as in Notsu and Tabata , Rui and Tabata , and Long and Yuan ).
Using the standard central difference along the numerical characteristic and the transformation given in Appendix, we have Therefore the BDF2 scheme takes the form and the CN scheme takes the form To compare the results, we need the first-step backward difference scheme (BDF1): Let denote the -error, where and .
Figure 1 illustrates convergence of BDF1 and BDF2 using numerical characteristics and using CI (Figure 1(a)) and CSI (Figure 1(b)), respectively. The straight lines are the first-order line and the second-order line, respectively. We can see that BDF2 exhibits higher-order convergence than BDF1.
Figure 2 shows the convergence of CN schemes using exact characteristics (Figure 2(a)) and using numerical characteristics (Figure 2(b)), respectively. We can see that when first-order numerical characteristics are involved, the computed result is poor. This confirms the fact that the higher-order numerical characteristic line is necessary to ensure reasonable convergence.
Figure 3 exhibits convergence of BDF2 with the smaller diffusion coefficient using numerical characteristics. We see that either using CI or CSI, convergence is approximately second-order.
Remark 5. As the particular cases of the general formulation, the following schemes are also tested: Both of them exhibit similar convergence properties as the BDF2 schemes.
In this paper we formulate a general SL scheme of second-order. By the MOL approach, this formulation includes all previously introduced schemes. We prove that this general scheme is second-order accurate even if the first-order characteristic line is involved. This formulation is very easy to extend to higher order cases. We see that MMOC can be considered a special version of SL method. Convergence properties of SL finite difference schemes are numerically tested.
Transformation of Derivatives
In this appendix we deal with the differentiation of composite functions. When a convection-diffusion equation is written to the Lagrangian form, directional derivatives appears. However, for spatial discretizations we need transformation of the directional derivatives.
Let us give the transformation of derivatives for the diffusion term resulting from SL method.
Let , and . By the chain rule, In the sequel, we simply write as . From (A.1) it follows that By differentiating both sides of (A.2) with respect to , and noting that we have Similarly, by differentiating both sides of (A.2) with respect to , and noting that we have Eliminating the mixed partial derivatives from the left-hand sides of both (A.5) and (A.7), we have where the omitted terms consist of the second derivatives of or with respect to , . If and are polynomials of of degrees not more than one (as the example in Section 5), these terms will become zeros.
Similarly, from (A.3), we have Manipulation of (A.8) and (A.9) leads to where the omitted terms consist of the second derivatives of or .
If are polynomials of of degrees not more than one, then from (A.10) it follows that Furthermore, if we use the exact characteristic line then from (A.11), it holds that If we use the numerical characteristic line then from (A.11), it holds that
This work is supported by the Major State Research Program of China Grant no. 1999032803, the Natural Science Foundation of Shandong Grants no. ZR2010AL021 and ZR2010AQ010, and the Project of Shandong Province Higher Educational Science and Technology Program Grant no. J11LA09. The authors would like to express their gratitude to Professor Yirang Yuan for the help and the reviewers for the valuable comments and opinions.
R. Courant, E. Isaacson, and M. Rees, “On the solution of nonlinear hyperbolic differential equations by finite differences,” Communications on Pure and Applied Mathematics, vol. 5, pp. 243–255, 1952.View at: Publisher Site | Google Scholar | Zentralblatt MATH | MathSciNet
A. Robert, “A stable numerical integration scheme for the primitive meteorological equations,” Atmosphere-Ocean, pp. 35–46, 1981.View at: Google Scholar
A. Staniforth and J. Côté, “Semi-Lagrangeian integration schemes for atmospheric models-A review,” Monthly Weather Review, vol. 119, pp. 2206–2223, 1991.View at: Google Scholar
P. K. Smolarkiewicz and J. A. Pudykiewicz, “A class of semi-Lagrangian approximations for fluids,” Journal of the Atmospheric Sciences, vol. 49, no. 22, pp. 2082–2096, 1992.View at: Google Scholar
J. Douglas Jr. and T. F. Russell, “Numerical methods for convection-dominated diffusion problems based on combining the method of characteristics with finite element or finite difference procedures,” SIAM Journal on Numerical Analysis, vol. 19, pp. 871–885, 1982.View at: Publisher Site | Google Scholar | Zentralblatt MATH | MathSciNet
C. N. Dawson, C. J. Van Duijn, and M. F. Wheeler, “Characteristic-galerkin methods for contaminant transport with nonequilibrium adsorption kinetics,” SIAM Journal on Numerical Analysis, vol. 31, no. 4, pp. 982–999, 1994.View at: Publisher Site | Google Scholar | Zentralblatt MATH | MathSciNet
J. Douglas Jr., R. E. Ewing, and M. F. Wheeler, “The approximation of the pressure by a mixed method in the simulation of immiscible displacement,” RAIRO Journal on Numerical Analysis, vol. 17, pp. 17–33, 1983.View at: Google Scholar | Zentralblatt MATH | MathSciNet
T. F. Russell, “Time stepping along characteristics with incomplete iteration for a Galerkin approximation of immiscible displacement in porous media,” SIAM Journal on Numerical Analysis, vol. 22, no. 2, pp. 970–1013, 1985.View at: Publisher Site | Google Scholar | Zentralblatt MATH | MathSciNet
R. Bermejo, “On the equivalence of semi-Lagrangian schemes and particle-in-cell finite element methods,” Monthly Weather Review, vol. 118, no. 4, pp. 979–987, 1990.View at: Google Scholar
E. Süli, “Convergence and nonlinear stability of the Lagrange-Galerkin method for the Navier-Stokes equations,” Numerische Mathematik, vol. 53, no. 4, pp. 459–483, 1988.View at: Publisher Site | Google Scholar | Zentralblatt MATH | MathSciNet
Y. Yuan, “The characteristic-mixed finite element method for enhanced oil recovery simulation and optimal order L2 error estimate,” Chinese Science Bulletin, vol. 38, pp. 1066–1070, 1993.View at: Google Scholar
R. E. Ewing and H. Wang, “A summary of numerical methods for time-dependent advection-dominated partial differential equations,” Journal of Computational and Applied Mathematics, vol. 128, no. 1-2, pp. 423–445, 2001.View at: Publisher Site | Google Scholar | Zentralblatt MATH | MathSciNet
G. Strang, “On the construction and comparison of different splitting schemes,” SIAM Journal on Numerical Analysis, vol. 53, pp. 506–517, 1968.View at: Publisher Site | Google Scholar | MathSciNet
R. E. Ewing and T. F. Russell, “Multistep Galerkin methods along characteristics for convection-diffusion problems,” in Advances in Computer Methods for Partial Differential Equations IV, R. Vichnevetsky and R. S. Stepleman, Eds., pp. 28–36, IMACS Rutgers University, New Brunswich, NJ, USA, 1981.View at: Google Scholar
H. Rui and M. Tabata, “A second order characteristic finite element scheme for convection-diffusion problems,” Numerische Mathematik, vol. 92, no. 1, pp. 161–177, 2002.View at: Publisher Site | Google Scholar | Zentralblatt MATH | MathSciNet
H. Notsu and M. Tabata, “A single-step characteristic-curve finite element scheme of second order in time for the incompressible Navier-Stokes equations,” Journal of Scientific Computing, vol. 38, no. 1, pp. 1–14, 2009.View at: Publisher Site | Google Scholar | Zentralblatt MATH | MathSciNet
K. Boukir, Y. Maday, and B. Métivet, “A high order characteristics method for the incompressible Navier-Stokes equations,” Computer Methods in Applied Mechanics and Engineering, vol. 116, no. 1–4, pp. 211–218, 1994.View at: Publisher Site | Google Scholar | Zentralblatt MATH | MathSciNet
K. Boukir, Y. Maday, B. Métivet, and E. Razafindrakoto, “A high-order characteristics/finite element method for the incompressible Navier-Stokes equations,” International Journal for Numerical Methods in Fluids, vol. 25, no. 12, pp. 1421–1454, 1997.View at: Publisher Site | Google Scholar | Zentralblatt MATH | MathSciNet
A. Bermúdez, M. R. Nogueiras, and C. Vázquez, “Numerical analysis of convection-diffusion-reaction problems with higher order characteristics/finite elements. Part I: time discretization,” SIAM Journal on Numerical Analysis, vol. 44, no. 5, pp. 1829–1853, 2006.View at: Publisher Site | Google Scholar | Zentralblatt MATH | MathSciNet
A. Bermúdez, M. R. Nogueirast, and C. Vázquez, “Numerical analysis of convection-diffusion-reaction problems with higher order characteristics/finite elements. Part II: fully discretized scheme and quadrature formulas,” SIAM Journal on Numerical Analysis, vol. 44, no. 5, pp. 1854–1876, 2006.View at: Publisher Site | Google Scholar | Zentralblatt MATH | MathSciNet
G. Fourestey and S. Piperno, “A second-order time-accurate ALE Lagrange-Galerkin method applied to wind engineering and control of bridge profiles,” Computer Methods in Applied Mechanics and Engineering, vol. 193, no. 39–41, pp. 4117–4137, 2004.View at: Publisher Site | Google Scholar
D. Xiu and G. E. Karniadakis, “A semi-Lagrangian high-order method for Navier-Stokes equations,” Journal of Computational Physics, vol. 172, no. 2, pp. 658–684, 2001.View at: Publisher Site | Google Scholar | Zentralblatt MATH | MathSciNet
D. Xiu, S. J. Sherwin, S. Dong, and G. E. Karniadakis, “Strong and auxiliary forms of the semi-Lagrangian method for incompressible flows,” Journal of Scientific Computing, vol. 25, no. 1, pp. 323–346, 2005.View at: Publisher Site | Google Scholar | Zentralblatt MATH | MathSciNet
M. Al-Lawatia, R. C. Sharpley, and H. Wang, “Second-order characteristic methods for advection-diffusion equations and comparison to other schemes,” Advances in Water Resources, vol. 22, no. 7, pp. 741–768, 1999.View at: Publisher Site | Google Scholar
M. Falcone and R. Ferretti, “Convergence analysis for a class of high-order semi-lagrangian advection schemes,” SIAM Journal on Numerical Analysis, vol. 35, no. 3, pp. 909–940, 1998.View at: Publisher Site | Google Scholar | Zentralblatt MATH | MathSciNet
R. Bermejo and J. Conde, “A conservative quasi-monotone semi-Lagrangian scheme,” Monthly Weather Review, vol. 130, no. 2, pp. 423–430, 2002.View at: Google Scholar
F. Xiao and T. Yabe, “Completely conservative and oscillationless semi-Lagrangian schemes for advection transportation,” Journal of Computational Physics, vol. 170, no. 2, pp. 498–522, 2001.View at: Publisher Site | Google Scholar | Zentralblatt MATH | MathSciNet
K. Toda, Y. Ogata, and T. Yabe, “Multi-dimensional conservative semi-Lagrangian method of characteristics CIP for the shallow water equations,” Journal of Computational Physics, vol. 228, no. 13, pp. 4917–4944, 2009.View at: Publisher Site | Google Scholar | Zentralblatt MATH | MathSciNet
E. Hairer and H. Wanner, Solving Ordinary Differential Equations I-Nonstiff Problems, Springer, 1987.View at: MathSciNet
J. D. Lambert, Computational Methods in Ordinary Differential Equations, John Wiley & Sons, London, UK, 1973.View at: MathSciNet
M. Zlamal, “Finite element methods for nonliear parabolic equations,” RAIRO Journal of Numerical Analysis, vol. 11, no. 1, pp. 93–107, 1977.View at: Google Scholar | Zentralblatt MATH | MathSciNet
X. Long and Y. Yuan, “Multistep characteristic method for incompressible flow in porous media,” Applied Mathematics and Computation, vol. 214, no. 1, pp. 259–270, 2009.View at: Publisher Site | Google Scholar | Zentralblatt MATH | MathSciNet