当前位置:首页 >> >>

Implicit-Explicit Methods for Time-Dependent PDE's


IMPLICIT-EXPLICIT METHODS FOR TIME-DEPENDENT PDE'S
URI M. ASCHER , STEVEN J. RUUTHy , AND BRIAN T.R. WETTONz

Abstract. Implicit-explicit (IMEX) schemes have been widely used, especially in conjunction with spectral methods, for the time integration of spatially discretized PDEs of di usion-convection type. Typically, an implicit scheme is used for the di usion term and an explicit scheme is used for the convection term. Reaction-di usion problems can also be approximated in this manner. In this work we systematically analyze the performance of such schemes, propose improved new schemes and pay particular attention to their relative performance in the context of fast multigrid algorithms and of aliasing reduction for spectral methods. For the prototype linear advection-di usion equation, a stability analysis for rst, second, third and fourth order multistep IMEX schemes is performed. Stable schemes permitting large time steps for a wide variety of problems and yielding appropriate decay of high frequency error modes are identi ed. Numerical experiments demonstrate that weak decay of high frequency modes can lead to extra iterations on the nest grid when using multigrid computations with nite di erence spatial discretization, and to aliasing when using spectral collocation for spatial discretization. When this behaviour occurs, use of weakly damping schemes such as the popular combination of Crank-Nicolson with second order Adams-Bashforth is discouraged and better alternatives are proposed. Our ndings are demonstrated on several examples.
region.

Key words. method of lines, nite di erences, spectral methods, aliasing, multigrid, stability AMS subject classi cations. 65J15,65M20

systems arising from spatially discretized time-dependent partial di erential equations. For problems with terms of di erent types, implicit-explicit (IMEX) schemes have been often used, especially in conjunction with spectral methods 7, 16]. For convectiondi usion problems, for example, one would use an explicit scheme for the convection term and an implicit scheme for the di usion term. Reaction-di usion problems can also be approximated in this manner. In this work we systematically analyze the performance of such schemes, propose improved new schemes and pay particular attention to their relative performance in the context of fast multigrid algorithms and of aliasing reduction for spectral methods. Consider a time-dependent PDE in which the spatial derivatives have been discretized by central nite di erences or by some spectral method. This gives rise to a large system of ODEs in time which typically has the form (1)

1. Introduction. Various methods have been proposed to integrate dynamical

u = f (u) + g(u) _

where kg k is normalized and is a nonnegative parameter. The term f (u) in (1) is some possibly nonlinear term which we do not want to integrate implicitly. This could be because the Jacobian of f (u) is non-symmetric and non-de nite and an
Department of Computer Science, University of British Columbia, Vancouver, BC, V6T 1Z4, Canada. The work of this author was partially supported under NSERC Canada Grant OGP0004306.
y Department of Mathematics, University of British Columbia, Vancouver, BC, V6T 1Z2, Canada. The work of this author was partially supported by an NSERC Postgraduate Scholarship. z Department of Mathematics, University of British Columbia, Vancouver, BC, V6T 1Z2, Canada. The work of this author was partially supported under NSERC Canada Grant OGP0122105. 1

1. INTRODUCTION

2

iterative solution of the implicit equations is desired, or the Jacobian could be dense, as in spectral methods, requiring the inversion of a full matrix at each time step. One may simply wish to integrate f (u) explicitly for ease of implementation. The term g (u), however, is a sti term which should be integrated implicitly to avoid excessively small time steps. Frequently g (u) is a linear di usion term, in which case the implicit equations form a linear system which is positive de nite, symmetric and sparse. Such systems can be solved e ciently by iterative techniques (e.g. 26, 11]). Thus, for problems of the form (1) it often makes sense to integrate g (u) implicitly and f (u) explicitly, yielding an IMEX scheme. The most popular IMEX scheme hitherto has been a combination of second order Adams-Bashforth for the explicit (\convection") term and Crank-Nicolson for the implicit (\di usion") term 16, 7, 1]. Applied to (1) this gives

un+1 ? un = 3 f (un ) ? 1 f (un?1 ) + g(un+1 ) + g(un)] k 2 2 2 where k is the constant discretization step size and un is the numerical approximation to u(kn). Other implicit-explicit methods have also been considered, e.g. 18, 25, 15].
(2)

A wide variety of other applications for IMEX schemes are also possible. For example, solutions to reaction-di usion systems arising in chemistry and mathematical biology can be computed using this technique. For these problems the nonlinear reaction term can be treated explicitly while the di usion term is treated implicitly. Examples of reaction-di usion systems from a biological standpoint can be found in 19]. We report about these in a separate paper 21]. Several authors have analyzed speci c IMEX schemes. For example, an experimental analysis of several implicit-explicit schemes including (2) was carried out in 2]. Some stability properties for certain second order IMEX schemes were determined in 25]. The schemes (17) and (24) below were considered in a Navier-Stokes context in 15]. These studies do not address how to choose the best IMEX scheme for a given system (1). In this paper, we seek e cient IMEX schemes for convection-di usion type problems using a systematic approach. First we derive, in x2, classes of linear multistep IMEX schemes. Classes of s-step schemes turn out to have optimal order s and to depend on s parameters. We then restrict our attention in x3 to a prototype advection-di usion equation in 1D

(3)

Ut = aUx + Uxx

where a and are constants, > 0, subject to periodic boundary conditions. A von-Neumann analysis (see, e.g. 23]) can then be applied, yielding in e ect a safe diagonalization of (1) and allowing us to consider a scalar test equation (4)

x = ( + i )x _

as is customary for ODEs (see, e.g. 12]), for values of , on an ellipse. An analysis of IMEX schemes for (4) is then performed, seeking methods which allow the largest stable time steps. In particular, for 1 we seek methods possessing a mild time step restriction since the system (1) is dominated by the implicitly handled \di usion

2. GENERAL LINEAR MULTISTEP IMEX SCHEMES

3

term". For other values of 0 we seek schemes which have reasonable time-stepping restrictions, e.g. comparable to the CFL condition 23]. In x4 we perform a variety of numerical experiments using nite di erence and pseudospectral methods on linear and nonlinear problems of convection-di usion type in one and two spatial dimensions. These experiments agree with the theory of x3. Furthermore, we discuss and demonstrate the use of IMEX schemes which yield strong decay of high frequency spatial modes. This property has important implications for the e ciency of time-dependent multigrid methods and of pseudospectral methods. An appropriate IMEX scheme (but not the popular (2)!) can reduce aliasing in pseudospectral methods. In the multigrid context, similar methods can reduce the number of multigrid cycles needed per time step, in e ect acting as smoothers (cf. 14, 10]). Conclusions and recommendations are summarized in x5. Perhaps the most surprising conclusion is that the most popular IMEX scheme (2) can essentially always be outperformed by other IMEX schemes. Moreover, a modi cation of (2) is almost always at least as good as (2) and at times much better. 2. General linear multistep IMEX schemes. We now derive s-step IMEX schemes for (1), s 1. Letting k be the discretization step size and un denote the approximate solution at tn = kn, these schemes may be written, (5)

k

X X 1 un+1 + 1 s?1 a un?j = s?1 b f (un?j ) +

k j=0

j

j =0

j

s?1 X j =?1

cj g(un?j )

where c?1 6= 0. See 8] for some stability and convergence results. For a smooth function u(t), expand (5) in a Taylor series about tn = n t to obtain the truncation error. This yields
X X 1 1 + s?1 a ]u(t ) + 1 ? s?1 ja ]u(t ) + (6) k j _ n j n
j =0 j =1

? ?

s?1 X

j =0 s?1 X

bj f (u(tn )) + k

s?1 X

j =1

jbj df jt=tn ? dt
s?1 X

? (p ? 1)!

j =?1 kp?1

cj g(u(tn)) ? k c?1 ?

j =?1 s?1 p?1 g X c?1 + cj (?j )p?1] d tp?1 jt=tn d j =1 s?1 X j =0

cj ] dg jt=tn ? dt

X kp?1 1 + s?1(?j )pa ]u(p)(t ) + p! j n j =1 p?1 p?1 s?1 X ? (pk? 1)! (?j )p?1bj d tp?f jt=tn d 1 j =1

+ O(kp ):

Applying (1) to the truncation error (6), an order p scheme is obtained provided 1+ 1?
s?1 X

aj

=
s?1 X j =0 s?1 X

0 =
s?1 X j =?1

s?1 X 1+ aj = ? jbj = c?1 ? jcj 2 j =1 2 j =1 j =1

j =1 s?1 j 2 X

jaj =

bj

cj

3. ANALYSIS FOR A TEST ADVECTION-DIFFUSION PROBLEM

4

(7)

X X X 1 + s?1 (?j )p a =s?1 (?j )p?1bj = c?1 + s?1 (?j )p?1cj : j p! j=1 p! (p ? 1)! j =1 (p ? 1)! j =1 (p ? 1)!

. . .

It is not di cult to prove the following theorem 22]. Theorem 2.1. For the s-step IMEX scheme (5), the following hold. (a) The 2p + 1 constraints of the system (7) are linearly independent, provided p s. Thus, there exist s-step IMEX schemes of order s. (b) An s-step IMEX scheme cannot have order of accuracy greater than s. (c) The family of s-step IMEX schemes of order s has s parameters We thus restrict all further discussion to s-step IMEX schemes of order s. 3. Analysis for a test advection-di usion problem. For the prototype problem (3), using the standard 2nd order centred approximations, D1 and D2, for the rst and second derivatives, respectively, we obtain the corresponding semi-discrete equations, _ Ui = aD1Ui + D2 Ui; 1 i M: (Here, a uniform spatial grid with M points has been employed.) Applying a discrete Fourier transform diagonalizes this system to xj = i j xj + j xj ; j = 1 : : :M _ where j and j are given by (8) ( ; ) = ( 2 cos(2 jh) ? 1]; a sin(2 jh))
j j

For notational convenience, we write (9) x = i x+ x _ We then consider values ( ; ) which lie on the ellipse of Figure 1. Note that the solution of (9) decays in time, jx(tn+1)j = e k jx(tn)j: Applying the general multistep IMEX scheme (5) to (9) yields (10)
X X X 1 xn+1 + 1 s?1 a xn?j = s?1 b i xn?j + s?1 c xn?j j j j k k
j =0 j =0 j =?1

h2

h

which is a linear di erence equation with constant coe cients. (In comparing (10) to (5) note that is buried in .) The solutions of the di erence equation (10) are of the form xn+1 = p1 1n + p2 2n + + ps sn where j is the j th root of the characteristic equation de ned by (11) (z ) (1 ? c?1 k)z s +
s?1 X j =0

(aj ? bj i k ? cj k)z s?j ?1

3. ANALYSIS FOR A TEST ADVECTION-DIFFUSION PROBLEM

5

β

a/h

?4ν/h2

Fig. 1

α . Half ellipse of possible ( ;

0

)

and pj is constant for j simple and a polynomial in n otherwise. Clearly, stability holds for j j j 1; j simple and j j j < 1 where j is a multiple root.

Remark. Below we will be interested in the entire range of

0. For the purely hyperbolic case = 0, our choice of the centred di erence D1 then leaves out the important characteristic-based methods. However, it is well-known that a one-sided scheme can be interpreted as a centred scheme plus an arti cial di usion term (e.g. 13]), so we do want to pay close attention to = O(h). 2 3.1. First order IMEX schemes. The one-parameter family of 1st order IMEX schemes for (1) can be written as (12) un+1 ? un = kf (un) + k (1 ? )g(un) + g(un+1)] and we restrict 0 1 to prevent large truncation error. The choice = 0 yields the Forward Euler scheme, un+1 ? un = kf (un ) + kg(un ): 1 This scheme is fully explicit, and will not be considered further. The choice = 2 yields the 2nd order Crank-Nicolson scheme when f = 0. Another possibility is to apply Backward Euler to g and Forward Euler to f . This choice ( = 1) yields, (13) un+1 ? un = kf (un) + kg(un+1 ): IMEX schemes such as (13) which apply a BDF discretization 9] to g and which extrapolate f to time step (n + 1) will be referred to as semi-implicit BDF (SBDF) schemes. Applied to the test equation (4), the general IMEX scheme (5) gives xn+1 = ( ; )xn where ( ; ) = 1 + k 1(1 ? ) + ik : ?k

3. ANALYSIS FOR A TEST ADVECTION-DIFFUSION PROBLEM

6

The stability region is thus f( ; ) : j ( ; )j 1g. Note that these one-step schemes trivially accommodate variable time-stepping and use relatively little storage. The choice = 1 is particularly attractive because strong decay occurs for large and negative. However, rst order IMEX schemes have the shortcoming that they are unstable near the nonzero -axis (i.e. for su ciently small and a 6= 0 xed but arbitrary) since

j (0; )j = j1 + ik j = 1 + k2 2 > 1:
Furthermore, at least a second order time-stepping scheme is often desirable since a second order spatial discretization is used. We thus consider higher order methods for the remainder of this paper. 3.2. Second order IMEX schemes. Approximating (1) to second order using IMEX schemes leaves two free parameters. If we centre our schemes in time about time step (n + ) to second order, we obtain the following family, 1 n+1 1 n?1 1 n n n?1 (14) k ( + 2 )u ? 2 u + ( ? 2 )u ] = ( + 1)f (u ) ? f (u )+ c c ( + 2 )g (un+1) + (1 ? ? c)g (un) + 2 g (un?1)]:

q

Some of these methods are quite familiar. For example, selecting ( ; c) = ( 1 ; 0) yields 2 the popular scheme (2), 1 un+1 ? un ] = 3 f (un ) ? 1 f (un?1 ) + g (un+1 ) + g (un)]: k 2 2 2 Because it applies Crank-Nicolson for the implicit part and second order AdamsBashforth for the explicit part, this scheme will be referred to as CNAB (CrankNicolson, Adams-Bashforth). Below, we show that the best asymptotic decay proper1 1 ties for = 2 occur when c = 8 . This choice gives 1 9 3 1 (15) k un+1 ? un ] = 3 f (un ) ? 1 f (un?1 ) + 16 g (un+1) + 8 g (un ) + 16 g (un?1 )] 2 2 Because of the obvious similarity to CNAB, this scheme will be called modi ed CNAB (MCNAB). Note that in comparison to CNAB, MCNAB does require the additional evaluation or storage of g (un?1). By setting ( ; c) = (0; 1) in (14) we obtain another method which has been applied to spectral applications (e.g. 4]), 1 un+1 ? un?1 ] = f (un ) + g (un+1) + g (un?1 )]: (16) 2k 2 This scheme uses leap frog explicitly and something similar to Crank-Nicolson implicitly (cf. 17]). For this reason, this method shall be referred to as CNLF (CrankNicolson, Leap Frog). Finally, setting ( ; c) = (1; 0) yields 1 n+1 n n?1 n n?1 n+1 (17) 2k 3u ? 4u + u ] = 2f (u ) ? f (u ) + g (u )

3. ANALYSIS FOR A TEST ADVECTION-DIFFUSION PROBLEM

7

It is easy to verify that at the origin, ( ; ) = (0; 0), the roots of are 1 and 2 ?1 . 2 +1 Thus, all of these schemes are stable at the origin, provided 0. Because the parameter c in (18) is always multiplied by we choose c according to stability properties for j j 1=k. For this case, the roots of the characteristic equation (18) are given approximately by which yields ( + c ) 2 + (1 ? ? c) + c = 0 2 2
p +c?1 (1 ? )2 ? 2c 1;2 = 2 +c

which shall be referred to as SBDF since this scheme is centred about time step (n +1). Other authors, e.g. 25], call this scheme extrapolated Gear. Having determined integration formulae, we direct our attention to obtaining stability properties for second order IMEX schemes. The second degree characteristic polynomial resulting from (14) applied to x = ( + i )x is given by _ (z ) = + 1 ? k( + c )]z 2 (18) 2 2 c 1 ? 2 + i k( + 1) + k(1 ? ? c)]z + ? 2 + i k ? k 2

For any ( ; c), evaluating (19)

D ;c max(j 1j; j 2j)

provides an estimate of the high frequency modal decay for large . Minimization over c determines the method with the strongest asymptotic high frequency decay for a particular . This yields (see 22]) 2 (20) c = (1 ? ) for > 0; 2 1 c 2 for = 0: The choice (21)

c= 1?

also proves useful (see 22]). The schemes SDBF, MCNAB and CNLF satisfy (20), the schemes SDBF and CNLF satisfy (21), but CNAB satis es neither. Also, minimization of (19) over and c determines that the SBDF scheme (17) possesses the strongest asymptotic decay of second order methods. Further information can be obtained from the stability contours in the ? plane. These plots are displayed in Figures 2 to 5. Figure 2 shows the contours for CNLF. This method is stable for all 0, provided k < h . Such a time step restriction is a undesirable since it applies even for large and small h. Furthermore, the decay of high frequency modes can be weak, tending to 1 as ! ?1. Comparison of CNLF to other second order methods such that = 0 suggests that CNLF produces the largest stability region among such methods.

3. ANALYSIS FOR A TEST ADVECTION-DIFFUSION PROBLEM
β

8
β

1/k

1/k

-10/k

Fig. 2. Stability contours for CNLF

α

0

-10/k

Fig. 3. Stability contours for CNAB

α

0

β

β

1/k

1/k

-10/k

Fig. 4. Stability contours for MCNAB

α

0

-10/k

Fig. 5. Stability contours for SBDF

α

0

The contours for CNAB are plotted in Figure 3. This method has a reasonable time step restriction for larger and small h. It is unstable near the -axis, however. It also su ers from poor decay of high frequency modes, since the decay tends to 1 1 as ! ?1. Using MCNAB, ( ; c) = ( 1 ; 1 ), the decay tends to 3 , a signi cant 2 8 improvement. See Figure 4 for these contours. The contours for SBDF are displayed in Figure 5. This method has the mildest time-stepping restriction when is large and h is small. The decay of high frequency modes is strong, tending to 0 as tends to ?1. This method, however, has the strictest time step limitation for small j j. We can now develop a quantitative method for describing time step restrictions for second order schemes. Such a method will help select which second order scheme to use for a particular problem. De ning the mesh Reynolds number 20] (22)

R ah

we plot in Figures 6 to 8 the theoretical time step restrictions for various and c. As can be seen from these gures, increasing allows larger stable time steps when R < 2. The case c = 1 ? also has the property that decreasing allows larger time steps for R > 2. Comparison of Figures 6 to 8 indicates that the largest time step can be applied using SBDF for R < 2 and CNLF for R > 2. This result physically corresponds to selecting SBDF when discrete di usion dominates, and CNLF when

3. ANALYSIS FOR A TEST ADVECTION-DIFFUSION PROBLEM

9

Fig. 6

. Time step restriction for various where c = 1 ?

=1/R

discrete convection dominates. From this perspective, the popular CNAB is only competitive when R 2. Frequently, an important consideration when choosing a second order scheme is what the constant of the truncation error is. For example, Crank-Nicolson is known to have a much smaller truncation error than second order BDF (see, e.g. 12]), so we expect CNAB to have a smaller truncation error than SBDF. The scheme MCNAB is expected to have a truncation error similar to CNAB, however. (Numerical experiments in x4 support this claim.) Because of its small truncation error and because it produces stronger high frequency spatial decay than CNAB, MCNAB may be preferred in certain problems over both CNAB and SBDF when R < 2. One disadvantage that all second order IMEX schemes with > 0 (i.e. essentially all except CNLF) share, is that their stability regions do not contain a portion of the -axis other than the origin. Speci cally, it was shown in 25] that when = 0 one of the roots, i , of (18) satis es

23]), providing no damping of high frequency error components anywhere. To obtain IMEX schemes which are stable for 0 and have strong decay for j j 1=k we must consider higher order schemes. 3.3. Third order IMEX schemes. Recall from x2 that we have a 3-parameter family of 3-step IMEX schemes of order 3. One possible parameterization yields, 1 1 1 1 (23) k ( 2 2 + + 3 + )un+1 + (? 3 2 ? 2 + 2 ? )un + 2

j i( k)j2 = 1 + ( 2 + 2 )( k)4 + > 1 for k su ciently small and > 0. Thus for all k, all second order schemes such that > 0 are unstable on the nonzero -axis. The CNLF scheme is well-known to be stable on the -axis provided k < h , but it is only marginally stable (see, e.g. a

3. ANALYSIS FOR A TEST ADVECTION-DIFFUSION PROBLEM

10

Fig. 7

2 . Time step restriction for various where c = (1?2 )

=1/R

Fig. 8

. Time step restriction for various where c = 0

=1/R

1 ( 3 2 + ? 1)un?1 + (? 1 2 + 6 )un?2 ] = 2 2 2+3 23 )f (un ) ? ( 2 + 2 + 4 )f (un?1 ) + ( 2 + + 5 )f (un?2 ) + ( 2 + 1 + 12 3 2 12 2+ 23 )g (un ) + ( 2 + c)g (un+1) + (1 ? 2 ? 3c + 12

3. ANALYSIS FOR A TEST ADVECTION-DIFFUSION PROBLEM

11

We now determine which methods produce the strongest asymptotic decay as ! ?1. For this case, the roots of the characteristic polynomial (25) are given approximately by 2+ 2 4 5 (26) ( 2 + c)z 3 + (1 ? 2 ? 3c + 23 )z 2 ? ( ? ? 3c + 3 )z + 12 ? c = 0: 12 2 By determining the solutions, 1; 2 and 3, of (26) we may evaluate D ; ;c max(j 1j; j 2j; j 3j)

2? 4 5 ( 2 + 3c ? 3 )g (un?1 ) + ( 12 ? c)g (un?2)]: These schemes are centred about time step (n + ) to third order, provided = 0. As for lower order schemes, the value of should be between 0 and 1 to avoid large truncation error. Also, the parameter c is multiplied by , so this parameter should be adjusted to modify large properties of a scheme. Setting ( ; ; c) = (1; 0; 0) yields 1 1 (24) k ( 11 un+1 ? 3un + 3 un?1 ? 3 un?2 ) = 3f (un ) ? 3f (un?1 ) + f (un?2 ) + g (un+1 ) 6 2 which is the third order BDF for the implicit part, and therefore is called third order SBDF. The third degree characteristic polynomial resulting from (23) applied to the test equation x = ( + i )x is given by _ 2+ 1 1 (25) (z ) = 2 2 + + 3 + ? ( 2 + c) k]z 3 2 23 23 3 ? 2 2 + 2 ? 1 + + ( + 3 + 1 + 12 )i k + (1 ? 2 ? 3c + 12 ) k]z2 2 2 3 2 + ? 1 + ( 2 + 2 + 4 )i k + ( ? 2 ? 3c + 4 ) k]z + 2 3 2 3 1 2 ? 1 + ( 2 + + 5 )i k + ( 5 ? c) k]: ? 2 6 2 12 12

to obtain an estimate of the high frequency model decay for large . Minimization over ( ; ; c) determines the method with the strongest asymptotic high frequency decay. Certainly if (27) D ; ;c = 0 then ( 0; 0; c0) minimizes the ampli cation as ! ?1. From (26) we satisfy (27) i 23 2 1 ? 0 ? 3c0 + 12 0 = 0; + +c 0 2 ? + 3c ? 4 0 3 0 = 0; 2 (28) + +c 0 2
0 0 0 2 0 0 2 0 0 2 0 0

5 12 0 ? c0 2+ 0 0 2 + c0

= 0:

3. ANALYSIS FOR A TEST ADVECTION-DIFFUSION PROBLEM
2 0 0 2 0 0

12

Each term is divided by ( + + c0 ) to allow the possibility of satisfying (27) by letting 2 + + c ) ! 1. Simpli cation of the system (28) yields ( 2 0 (3 0 ? 1)( 0 ? 1) = 0; (29) + +c 0 2 2 ? 6( 0 ? 0) = 0; 0 (30) + +c 0 2
2 0 0 2 0 0

(31)
2 0 0

For ( + + c0) nite there are two possibilities, ( 0; 0; c0) = (1; 0; 0) and ( 0; 0 ; c0) = 2 1 ; ? 4 ; ? 5 ), both of which specify third order SBDF. For ( + + c0 ) in nite, (29) (3 3 9 2 implies j 0j jc0j. Using this in (30) implies j 0 j jc0j. Applying these results to (31) results in a contradiction, so ( + + c0) must be nite. 2 Because third order SBDF has the strongest asymptotic decay of third order IMEX schemes, special attention is given to its properties throughout the remainder of this section. All schemes considered here are stable on a segment of the -axis including the origin ( ; ) = (0; 0), as can be veri ed by an analysis of (25) 22]. Further information about stability in the - plane can be obtained by plotting maxfjz j : (z ) = 0g, where (z ) is de ned in (25). These stability contours are displayed in Figures 9 to 14. We begin by examining if it is possible to arrive at a stable scheme for any xed . For a xed, but arbitrary and for j j; jcj ! 1 we obtain an approximate local minimum of D ; ;c if c = 0:4661: Using these parameters, the scheme simpli es to 1 5 (32) k (un+1 ? un ) = 23 f (un ) ? 4 f (un?1 ) + 12 f (un?2 ) + :4661 g (un+1) 12 3 +:5184 g (un) + :0650 g (un?1) ? :0494 g (un?2) which applies third order Adams-Bashforth to the explicit term. (Note that disappears in (32) through the limiting process.) The stability contours of Figure 9 suggest that this method is stable for all 0 provided k < 0:62 h . This restriction is more a severe than that for third order Adams-Bashforth applied to x = i x because of the _ dip in the stability contours when < 0. As mentioned previously, the -axis is included in the absolute stability region. This result demonstrates that for third order methods, it is possible to nd methods for any which are stable for all 0 by varying and c. In x3.2, the most interesting second order schemes were produced by selecting equal to 0, 1 or 1. We consider schemes for these values of below. The parameters 2 and c are chosen to produce schemes which allow large stable time steps as ! 1. For example, the method ( ; ; c) = (0; ?2:036; ?:876) of Figure 10 is stable for all 0 provided k 0:67 h . Similarly ( ; ; c) = (:5; ?1:21; ?:5) of Figure 12 is stable a for all 0 if k 0:65 h . In both these cases substantially larger time steps can be a taken for large or moderate j j than for the method (32). Furthermore, stronger high frequency decay occurs for these methods than for method (32). Recall that the strongest decay as ! ?1 occurs for third order SBDF. The stability contours for this method are shown in Figure 13. This plot together with the
2 0 0 2 0 0

5 12 0 ? c0 2 0+ 0 2 + c0

= 0:

3. ANALYSIS FOR A TEST ADVECTION-DIFFUSION PROBLEM
β

13
β

1/k

1/k

-9/k

Fig. 9. Stability contours for ( ; ; c) = lim ( ; ; :4661 ) !1

α

0

-9/k

Fig. 10. Stability contours for ( ; ; c) = (0; ?2:036; ?:876)
β β

α

0

1/k

1/k

-9/k

Fig. 11. Stability contours for ( ; ; c) = ( 1 ; 0; 0) 2
β

α

0

-9/k

Fig. 12. Stability contours for 1 ( ; ; c) = ( 2 ; ?1:21; ?:5)
β

α

0

1/k

1/k

-9/k

Fig. 13. Stability contours for 3rd order SBDF, ( ; ; c) = (1; 0; 0)

α

0

-.2/k

Fig. 14. Zoom-in around = 0 for 3rd Order SBDF

α

0

0 .1/k

zoom-in of Figure 14 suggest that third order SBDF is stable for all 0 provided h . The plot of Figure 13 also indicates that very large time steps can be taken k < 0:62 a for large or moderate j j. Although applying = 1 and nonzero and c may allow somewhat larger stable time steps we focus on = c = 0 since other choices degrade the strong asymptotic decay. For 6 1, the choice = c = 0 is not recommended. As seen in Figure 11, 1 ( ; ; c) = ( 2 ; 0; 0) results in a small stability region. This plot indicates that very small time steps are needed for moderate or large h , since the ellipse of Figure 1 must lie in the stable region.
2

4. FURTHER CONSIDERATIONS AND NUMERICAL EXPERIMENTS

14

4th order, 4-step IMEX scheme is a 4 parameter family of methods. From the previous section we anticipate that the 4th order SBDF may have good stability properties. For u = f (u) + g(u), this scheme is given by _ 1 ( 25 un+1 ? 4un + 3un?1 ? 4 un?2 + 1 un?3 ) = (33) k 12 3 4 4f (un ) ? 6f (un?1 ) + 4f (un?2 ) ? f (un?3 ) + g (un+1 ): The characteristic polynomial obtained from applying (33) to x = ( + i )x is _ 25 4 (z ) = ( 12 ? k)z 4 ? (4 + 4i k)z 3 + (3 + 6i k)z 2 ? ( 3 + 4i k)z + 1 + i k: 4 Contour plots similar to those in Figures 13 and 14 were obtained for the 4th order SBDF as well (see 22] for these plots). The scheme is stable for k 0:52 h and larger a time steps are permitted as increases. However, the stability region is smaller than for the third order case, so smaller time steps may be required. Furthermore, the axis is closer to the boundary of the stability region for fourth order SBDF, suggesting that third order SBDF may dissipate error better for 0. Third and fourth order SBDF methods are good choices for IMEX schemes for some problems. For all 0, these methods are stable for a time step restriction similar to the CFL condition. Greater accuracy and strong high frequency decay also make these methods very attractive. Nonetheless, for many problems second order methods are preferred. Higher order methods require more storage, and more work per time step. This extra expense could be justi ed if greater accuracy permitted larger time steps. Third and fourth order schemes, however, have more severe time step restrictions than second order schemes. In fact, Figure 15 shows that larger stable time steps can be taken with second order SBDF when R < 9 for the linear advection-di usion problem using second order spatial discretizations. CNLF allows larger stable time steps than either third or fourth order SBDF for R > 1. Third order SBDF should be e cient for problems which require strong decay for j j 1=k and a moderate time step restriction for R > 10. It should also be e ective for problems where R 1, since a portion of the -axis is within the stability region. Fourth order SBDF has a particularly severe time step restriction for the advectiondi usion problem when R is moderate or small. For example, when R = :125, fourth order SBDF can only apply one tenth the time step of second order SBDF, as can be seen from Figure 15. This restriction on the step size would appear to limit fourth order SBDF to problems where accuracy, and not stability, is the reason for limiting the time step size. 4. Further considerations and numerical experiments. The previous section has dealt with stability properties of IMEX schemes for the one-dimensional linear constant coe cient advection-di usion equation. These results provide necessary, but not su cient conditions for stability for variable coe cient and nonlinear convectiondi usion problems. In this section we summarize numerical experiments which verify that our analysis for the simple, linear problem can be useful for determining which IMEX scheme to apply to more complicated problems. In addition, strongly damping schemes are shown to be more e cient in certain spectral collocation and multigrid applications.

3.4. Fourth order IMEX schemes and a comparison. From x2, the general

4. FURTHER CONSIDERATIONS AND NUMERICAL EXPERIMENTS

15

Fig. 15

. Time step restrictions for various IMEX schemes

=1/R

4.1. Finite di erence approximations in 1D. Example 1. To examine nonzero viscosity behaviour, we consider the one-dimensional

In order to calculate starting values for multistep IMEX schemes, we use one-step (low order) IMEX schemes with a very small time step, unless otherwise noted.

variable coe cient problem (34) ut + sin(2 x)ux = uxx subject to periodic boundary conditions on the interval 0,1] and initial condition u(x; 0) = sin(2 x): We use centred second order di erences for the spatial derivatives. Note that the solution is smooth for all 0: To test the theory's predictions for small mesh Reynolds numbers (22), the model 1 problem (34) was approximately solved using discretization step sizes h = 63 in space and k = 1:8h in time. Use of such step sizes is appropriate only for strongly damped ow. Utilizing several IMEX schemes, computations to time t = 2 are performed for viscosities, , in the range :01 :1. These values correspond to mesh Reynolds numbers, R, in the range 1:59 R 0:159. From Figures 6, 7, 8 and 15 the theory predicts that for these step sizes SBDF is stable for a larger viscosity interval than any other scheme. As decreases, third order SBDF followed by CNAB should become unstable. Stability for MCNAB should be similar to, and somewhat better than, CNAB. Furthermore, fourth order SBDF and CNLF should be unstable over the entire interval because the theoretical stability restriction is violated. 1 By comparing results to those for h = 504 and k = :225h using SBDF, the max norm relative errors for second order schemes are evaluated. Third and fourth order

4. FURTHER CONSIDERATIONS AND NUMERICAL EXPERIMENTS

16

Fig. 16

. Large viscosity behaviour for various IMEX methods

schemes are compared to computations with these same step sizes but using third and fourth order SBDF. The resultant errors are plotted against in Figure 16. Fourth order SBDF and CNLF are not included because they are unstable over the entire interval. The plots of the gure clearly coincide with the results of the theory. Figure 16 also indicates that SBDF is the only stable method for the above choice of discretization parameters when 0:0015 < < 0:0025. (The values of where the curves turn upwards correspond to the onset of instability for the given values of h and k.) This agrees with the prediction that SBDF allows the largest stable time steps for small mesh Reynolds numbers (see x3.2). Although third order SBDF has a smaller stability interval, it may be useful in problems where high accuracy is needed since it produces a smaller error than second order methods when stable. To test the theory's predictions for large mesh Reynolds numbers, example (34) 1 was solved using discretization step sizes h = 81 in space and k = :9h in time. Using several IMEX schemes, computations to time t = 2 are performed for viscosities, , in the range :001 :01. These values correspond to mesh Reynolds numbers, R, in the range 12:3 R 1:23. From Figures 6, 7, 8 and 15 the theory predicts that for these step sizes only CNLF is stable over the entire viscosity interval. As decreases, third order SBDF followed by SBDF and nally CNAB should become unstable. Stability for MCNAB should be similar to CNAB. Furthermore, fourth order SBDF should be unstable over the entire interval because the theoretical stability restriction is violated when < :05 for these step sizes. 1 By comparing results to those for h = 648 and k = :225h using SBDF, the max norm relative errors for second order schemes are evaluated for these computations. Third and fourth order schemes are compared to computations with these same step sizes but using third and fourth order SBDF. The resultant errors are plotted against

4. FURTHER CONSIDERATIONS AND NUMERICAL EXPERIMENTS

17

Fig. 17

. Small viscosity example for various IMEX methods

in Figure 17. Fourth order SBDF is not plotted because it is unstable over the entire interval. The plots of the gure support the results of the theory. Figure 17 also indicates that CNLF is the only stable method for < 0:002. This agrees with the prediction that CNLF allows the largest stable time steps for large mesh Reynolds numbers (see Section 3.2). CNLF is a particularly attractive choice because it has the smallest error of second order methods. Use of SBDF is not recommended because it has the smallest stable interval and the largest error among any of the second order scheme considered. For the same problem using k = :5h, dramatically di erent results are predicted because all the methods satisfy their theoretical stability restrictions. Indeed, computations for CNLF, third and fourth order SBDF all produce errors which nearly coincide, since spatial error dominates the solutions. Other second order methods appear stable and produce only slightly less accurate results. 2 problem (35)

Example 2. To examine the limit = 0, we consider the one-dimensional nonlinear
1 ut + 2 cos(2 t)(1 + u)ux = 0

having periodic boundary conditions on the interval 0,1] and initial conditions

u(x; 0) = sin(2 x):
As in the previous two sections, we use second order centred di erences to ap1 proximate ux . For h = 80 and k = :5h, computations are performed to time t = 100. The time step value k = :5h was used to ensure that fourth order SBDF satis ed the

4. FURTHER CONSIDERATIONS AND NUMERICAL EXPERIMENTS

18

Fig. 18

. Burgers equation for various t for = 0:01

stability restriction k :52h. Using the fact that the exact solution to this problem at integer t equals the initial data, i.e. u(x; n) = sin(2 x); n = 0; 1; 2; : : : we compute the error in the solution at time t = n; n = 1; 2; : : :; 100. The results again agreed with the predictions from x3: all second order schemes tested, such that > 0, were unstable, while CNLF, third order SBDF and fourth order SBDF were all stable. The SBDF schemes are dissipative as well. For fuller details, see 22]. 2 (36) ut + uux = uxx subject to periodic boundary conditions on the interval -1,1] and initial conditions u(x; 0) = sin( x): A plot of the solution of the Burgers equation for = :01 at several di erent times is provided in Figure 18. This computation was produced using Chebyshev collocation with 40 basis functions and k = 1=160 using SBDF. The next few subsections discuss Fourier and Chebyshev collocation implementations for the above model problem. See 7] or 3] for details on these methods. 4.2.1. Fourier spectral collocation. Since the problem of Example 3 is periodic, we expect that Fourier series makes a good basis of trial functions for this problem. Indeed, since 2 2 ut = ? 1 @u + @ u 2 @x @x2

4.2. Spectral collocation. Example 3. For our next application we consider the well-known Burgers equation,

4. FURTHER CONSIDERATIONS AND NUMERICAL EXPERIMENTS

19

we see that ut is antisymmetric for u antisymmetric. By selecting an initial condition which is antisymmetric we guarantee that u remains antisymmetric for all t. Since only these components of the series contribute to the solution we use a Fourier sine series. We thus approximate u by the series

uN (x; t) =

N X j =1

j (t)sin(jx):

Using initial conditions 1 (0) = 1 and k (0) = 0; k 6= 1, the j (t) are determined by j enforcing the di erential equation at the collocation points, xj = N ? 21 ; 1 j N , N i.e.
"

@uN + u @uN = @ 2uN N @t @x @x2

#

x=xj

:

This scheme is also called a pseudospectral method since the nonlinear convection term is evaluated in physical space. By applying Fourier sine collocation with 40 basis functions and k = 1=40, we solve the model problem at each time step to time t = 2. Because the system is small, the implicit equations are solved in physical space using LU decomposition. In larger systems where greater e ciency is needed these would be solved in Fourier space using transform methods (see 7]). For CNAB, MCNAB, SBDF and third order SBDF the max norm relative error is plotted against viscosity (see Figure 19). The 1 \exact solution" is based on a computation using N = 80 modes and k = 3200 and third order SBDF. CNLF was not included because the theoretical stability restriction is violated. This can be easily seen because the linear advection-di usion equation has eigenvalues (ian ? n2 2) for eigenfunctions ein x . From these eigenvalues we know that the stability restriction is k < a 1n , which is violated initially because ju(x; 0)j1 = 1. As expected, third order SBDF has the smallest error of any of the methods when stable. The stable region, however, is smaller than that for CNAB or SBDF. CNAB and MCNAB are once again very similar in behaviour with the modi ed version being marginally better. SBDF appears to allow the largest stable time step when 0:01. 1 Further re nement of the time step to k = 160 leaves third order SBDF as the method of choice over the entire interval. Such a re nement may be unnecessary in this example because the error is less than 1% for a step size which is 4 times larger. 4.2.2. Chebyshev spectral collocation. Because the solution to the problem is periodic and anti-symmetric we know that u( 1; t) = 0 for all t. Using this fact, we solve (36) subject to the homogeneous Dirichlet boundary conditions u( 1; t) = 0, using a pseudospectral Chebyshev collocation scheme. The Gauss-Chebyshev grid,

xi = cos (2i ? 1) =2N ]; i = 1; : : :; N
is used for collocation points. Similar to the Fourier case, the max norm relative errors are evaluated for several IMEX schemes using k = 1=40 and N = 40. As can be seen from Figure 20, the results are qualitatively similar to those of the Fourier case. SBDF performs particularly well for smaller viscosities. From the theory, it has the widest stability region for large j j

4. FURTHER CONSIDERATIONS AND NUMERICAL EXPERIMENTS

20

Fig. 19

. Fourier spectral collocation for Burgers equation

Fig. 20

. Chebyshev spectral collocation for Burgers equation

and thus is best able to accommodate the rapidly growing eigenvalues associated with Chebyshev collocation. Both Chebyshev and Fourier collocation can be a ected by aliasing 7]. We consider aliasing for the Fourier case next.

4. FURTHER CONSIDERATIONS AND NUMERICAL EXPERIMENTS

21

terms produce frequencies that cannot be represented in the basis, and thus contribute erroneously to lower frequencies. For instance, in Example 3, the Fourier sine mode sin(mx), when explicitly evaluated in the convection term produces the contribution sin(mx) @ sin(mx) = m sin(2mx): @x 2 If 2m is greater than the number of Fourier sine basis functions, N , this frequency cannot be represented correctly and aliasing occurs. We now proceed to demonstrate that this behaviour can plague poorly spatially resolved computations when applying weakly damping IMEX schemes. We compute solutions for the model Burgers equation (36) subject to periodic boundary conditions and the initial conditions of Example 3 modi ed to read

4.2.3. Aliasing in pseudospectral methods. Aliasing occurs when nonlinear

u(x; 0) = 0:98 sin(2 x) + HF (x)
where

HF (x) 0:01 sin(61 x) + 0:01 sin(62 x):
We use Fourier sine collocation with N = 64 basis functions, and integrate to time t = 2 with a viscosity = 0:1. The value of the approximate solution at time t = k is obtained using rst order SBDF with the same step size. To represent the type of high frequency information that could be produced during a computation, we add HF (x) to the solution at t = k. This ensures that high frequency information remains after the strongly damping rst order SBDF step. Subsequent steps are then computed using the relevant second order scheme. For third order SBDF, the value at t = 2k is also needed. For the purpose of demonstrating aliasing e ects, this value is computed using CNAB since we wish to retain most of the added high frequencies. The max norm relative errors for several IMEX schemes are computed by comparing results to those for SBDF using N = 128 and k = 241N : In the third order case, results are compared to those for third order SBDF using N = 128 and k = 241N : These errors are summarized below. CNLF CNAB MCNAB SBDF 3rd order SBDF 1 .575 .0056 .0072 .0045 64 .92 1 .060 .022 .0013 .00079 .0010 192 1 For the case k = 64 we note that the error for CNAB is far greater than for MCNAB, SBDF or third order SBDF. Using CNAB, a non-aliased computation using 1 the 2/3's rule 7] and k = 64 results in a relative error of less that 10?3. This nonaliased result, along with its aliased counterpart, are plotted in Figure 21. The \error" curve in this gure is that of the aliased CNAB. From the gure it is clear that the main component of the error is proportional to sin( x). This low frequency mode is not part of the exact solution, and must be an aliasing artifact from high frequency components. It is interesting to note that even after 128 time steps the numerical solution is still plagued by high frequency components which have not yet decayed. A similar study of CNLF reveals that it su ers from aliasing as well.

k

4. FURTHER CONSIDERATIONS AND NUMERICAL EXPERIMENTS

22

Fig. 21

. Aliased and non-aliased computations for CNAB

1 Resorting to a smaller time step, k = 192 , makes a substantial improvement in the solution for CNAB and CNLF. Even so, the aliasing error for CNLF is su ciently large that further re nement is likely required. We conclude that use of a strongly damping scheme such as SBDF, MCNAB or third order SBDF gives an inexpensive method to reduce aliasing in poorly resolved computations. Application of weakly damping schemes like CNAB and especially CNLF may necessitate undesirably small time steps to produce the high frequency decay needed to prevent aliasing in an aliased computation. Alternatively, weakly damping schemes may be used in conjunction with an anti-aliasing technique such as the 3/2's rule or the 2/3's rule 7]. However, these anti-aliasing techniques either increase the expense of the computation or produce a severe loss of high frequency information at each time step. 2

4.3. Time-dependent multigrid in two spatial dimensions. Example 4. The next model that we consider is the 2D convection-di usion problem, (37) ut + (u r)u = u where u (u; v ). This model incorporates some of the major ingredients of the 2D
Navier-Stokes equations. Indeed, it is often being solved as part of projection schemes for incompressible ows (see 13]). We carry out our computations on the square

the numerical solution is arti cial. However, this allows us to study the propagation of errors related to these e ects in simple examples within a controlled environment. In general, high frequency solution and error components are generated at each time step by nonlinear terms, forcing terms and boundary e ects (cf. 6]). 2

Remark. It may be argued that our addition of the high frequency components to

4. FURTHER CONSIDERATIONS AND NUMERICAL EXPERIMENTS

23

0,1] 0,1] and consider periodic boundary conditions and the initial conditions

u(x; y; 0) = sin 2 (x + y)]; v(x; y; 0) = sin 2 (x + y)]:
For spatial discretization we use standard second order centred di erences, as before. For IMEX schemes the convection term, (u r)u, is handled explicitly and the di usion term, u, is handled implicitly. This treatment yields a positive de nite, symmetric, sparse linear system to solve at each time step. Such systems are solved e ciently using a multigrid algorithm, the components of which are outlined next. To solve the implicit equations we apply a Full Multigrid (FMG), Full Approximation Scheme (FAS) algorithm 5]. Rather than applying the standard algorithm at each time step, we use the ideas of 6]. Smoothing is accomplished using Red-Black Gauss-Seidel. This relaxation technique is chosen because it has a very good smoothing rate. Prolongation is accomplished using bilinear interpolation, and restriction by full weighting. The standard FMG cycle is modi ed so that the rst coarse grid correction is performed before any ne grid relaxation. Based on the assumption that the increment between time steps is smooth, ne grid relaxation is most e ective after the smoothness of the increment is accounted for (see 6]). Interpolation for the FMG step is bilinear. In 6] it is argued that for time-dependent problems higher order interpolation only decreases high frequency error, which tends to dissipate in parabolic systems anyway. This suggestion is utilized here; indeed experiments with cubic interpolation did not produce any reduction in the number of multigrid iterations to achieve a given tolerance. Another recommendation of 6] is to avoid the nal smoothing pass at each time step, in order to reduce aliasing from Red-Black Gauss-Seidel relaxation. Aliasing is not a major source of error in Example 4, however, so this suggestion was not utilized. For this problem we use a residual test with a tolerance TOL to determine the number of ne grid iterations to perform at each time step. Next we show that the choice of time-stepping scheme can a ect the number of ne grid iterations to achieve a given residual tolerance. 1 The model problem (37) was approximated using step sizes h = 128 and k = 0:00625 and residual tolerance TOL=0.003. After the rst time step, high frequency information

HF (x; y) = 0:005 cos 2 (64x + 63y)]
was added to each of u and v , to represent the type of high frequency information that might be produced during a computation. For several second order IMEX schemes, the average number of ne grid iterations at each time step was computed. The result using V-cycles is graphed in Figure 22. Strongly damping schemes such as SBDF and MCNAB require roughly 1 iteration per time step. Weakly damping schemes required far more e ort to solve the implicit equations accurately, because lingering high frequency components necessitate more work on the nest grid. This is evident since CNAB requires more than 2 iterations per time step and CNLF required 3. Using a W-cycle improves the relative e ciencies of CNLF and CNAB. Even for these cases, however, nearly twice the number of ne grid iterations were required

5. CONCLUSIONS AND RECOMMENDATIONS

24

3

Number of Iterations per Time Step

2 CNLF 1 CNAB MCNAB 0
1 1 (0,1) ( 2 ; 0) (1; 8) 2 Second Order Method ( ; c)
Fig. 22

SBDF (1,0)

1 . Multigrid V-cycle iterations, = 0:03; h = 128

than for more strongly damping schemes such as SBDF and MCNAB. See 22] for a plot. 1 Even for a smaller viscosity, = 0:02, and a much coarser mesh h = 32 , the performance of CNLF su ers. It uses about 30% more iterations than SBDF or MCNAB to achieve the desired tolerance. This result uses TOL=.009 and is plotted in Figure 23. Thus we conclude that use of a poorly damping IMEX scheme such as CNAB or especially CNLF can necessitate extra iterations on the nest grid in multigrid solutions to the implicit equations for small mesh Reynolds numbers, R < 2. For large mesh Reynolds numbers, R > 2, this e ect was not observed. 2 5. Conclusions and recommendations. IMEX schemes are not a universal cure for all problems. It is not di cult to imagine situations where a fully implicit or fully explicit scheme is preferable. However, these schemes can be very e ective in many situations, some of which are depicted in this paper. It may be important to choose an IMEX scheme carefully. The usual CNAB can be signi cantly outperformed by other IMEX schemes. Based on observations in this and previous sections we provide a few guidelines for selecting IMEX schemes for convection-di usion problems. 5.1. Finite di erences. For the nite di erence ( nite element, nite volume) case, begin by determining an estimate for the mesh Reynolds number, R = ah , where represents viscosity, a represents characteristic speed and h grid spacing. For problems where R 2 application of CNLF, or a third or fourth order scheme is reasonable. Third and fourth order SBDF can be applied to these problems, since a portion of the -axis is included in the stability region even though these methods were selected primarily for their large properties. Of the methods we have considered in detail, CNLF has the mildest time step restriction and is non-dissipative while third

5. CONCLUSIONS AND RECOMMENDATIONS

25

2

Number of Iterations per Time Step

1

CNLF CNAB MCNAB SBDF

0

1 1 (0,1) ( 2 ; 0) (1; 8) 2 Second Order Method ( ; c)
Fig. 23

(1,0)
1 32

. Multigrid W-cycle iterations, = 0:02; h =

order SBDF is the most dissipative. All other second order schemes should be avoided in this case. Explicit schemes should also be considered for these problems. For R > 2 use of CNLF or third order SBDF appears appropriate. CNLF has the mildest time step restriction, but accuracy concerns could make third order SBDF competitive. Avoid use of SBDF in this case. For problems of this type, a study to determine when explicit schemes are competitive would be interesting. For R 2 the theory predicts that many second order IMEX schemes have similar time step restrictions. A study to determine the method with the smallest truncation error would be useful in this case. For greater accuracy, third order SBDF appears to be more useful than fourth order SBDF, since its time step restriction is less severe. If strong decay of high frequency spatial modes is a desirable characteristic then CNLF should be avoided. For R < 2, use of SBDF permits the largest stable time steps. The modi ed CNAB scheme, MCNAB, can also be applied to problems of this type, although its time step restriction is somewhat stricter. Third order SBDF is recommended when high accuracy is needed. Numerical experiments in Section 4.3 demonstrate that in multigrid solutions to the implicit equations, application of a strongly damping method is prudent. MCNAB or SBDF should be useful in such problems. Avoid use of CNLF in this case. 5.2. Spectral methods. Because the eigenvalues for spectral methods are different than for nite di erences as is their meaning (see 24]), we cannot expect the stability time step restrictions from x3 to hold quantitatively. For this reason, an analysis of the linear advection-di usion equation for Chebyshev and Fourier spectral methods would be interesting. However, this is outside the scope of this paper. Numerical experiments for the Burgers equation were made for small to moderate mesh Reynolds numbers. For these problems, CNLF has a very severe time step

REFERENCES

26

restriction. These computations also suggest that SBDF has the mildest time step restriction of the methods considered. This result is particularly pronounced in the Chebyshev collocation case. The large mesh Reynolds number case for spectral collocation was not considered. In this case, a comparison of the relative e ciencies of IMEX schemes and fully explicit schemes would be interesting. Third order SBDF appears to be an e cient method for problems where high accuracy is needed. In problems where aliasing occurs, a strongly damping scheme, such as SBDF, MCNAB or third order SBDF can be used to inexpensively reduce aliasing. Application of a weakly damping scheme such as CNAB or CNLF in poorly spatially resolved, aliased computations should be avoided.
REFERENCES 1] L. Abia and J.M. Sanz-Serna. The spectral accuracy of a fully-discrete scheme for a nonlinear third order equation. Computing, 44:187{196, 1990. 2] C. Basdevant, M. Deville, P. Haldenwang, J.M. Lacroix, J. Ouazzani, R. Peyret, P. Orlandi, and A.T. Patera. Spectral and nite di erence solutions of the Burgers equation. Computational Fluids, 14:23{41, 1986. 3] J.P. Boyd. Chebyshev & Fourier Spectral Methods. Springer-Verlag, 1989. 4] M.E. Brachet, D.I. Meiron, S.A. Orszag, B.G. Nickel, R.H. Morf, and U. Frisch. Small-scale structure of the Taylor-Green vortex. J. Fluid Mechanics, 130:411{452, 1983. 5] A. Brandt. Guide to multigrid development. In W. Hackbusch and U. Trottenberg, editors, Multigrid Methods, pages 220{312. Springer-Verlag, Berlin, 1982. 6] A. Brandt and J. Greenwald. Parabolic multigrid revisited. In Proc. 3rd European Conference on Multigrid Methods, Bonn, Oct. 1990. 7] C. Canuto, M.Y. Hussaini, A. Quarteroni, and T.A. Zang. Spectral Methods in Fluid Dynamics. Springer-Verlag, 1987. 8] M. Crouzeix. Une methode multipas implicite-explicite pour l'approximation des equations d'evolution paraboliques. Numer. Math., 35:257{276, 1980. 9] C.W. Gear. The automatic integration of sti ordinary di erential equations. In Proc. IFIP Congress, pages A81{85, Ljubjana, Yugoslavia, 1986. 10] B. Gustafsson and P. Lotstedt. Analysis of the multigrid method applied to hyperbolic equations. Manuscript, Dept. Scient. Comp., Uppsala, 1992. 11] W. Hackbusch. Multi-Grid methods and Applications. Springer-Verlag, 1985. 12] E. Hairer, S.P. Norsett, and G. Wanner. Solving Ordinary Di erential Equations I. SpringerVerlag, 1987. 13] C. Hirsch. Numerical Computation of Internal and External Flows. Wiley, 1988. 14] A. Jameson. Computational transonics. Comm. Pure Appl. Math., XLI:507{549, 1988. 15] G.E. Karniadakis, M. Israeli, and S.A. Orszag. High-order splitting methods for the incompressible Navier-Stokes equations. J. Computational Phys., 97:414{443, 1991. 16] J. Kim and P. Moin. Application of a fractional-step method to incompressible Navier-Stokes equations. J. Computational Phys., 59:308{323, 1985. 17] H.O. Kreiss. Numerical methods for solving time-dependent problems for partial di erential equations. Lecture notes, Univ. de Montreal, 1978. 18] P.C. Meek and J. Norbury. Two-stage, two-level nite di erence schemes for nonlinear parabolic equations. IMA J. Num. Anal., 2:335{356, 1982. 19] J.D. Murray. Mathematical Biology. Springer-Verlag, 1989. 20] R. Peyret and T. Taylor. Computational Methods for Fluid Flow. Springer-Verlag, 1983. 21] S. Ruuth. Implicit-explicit methods for reaction{di usion problems. In preparation, 1993. 22] S. Ruuth. Implicit-explicit methods for time-dependent PDE's. Master's thesis, Ins. Appl. Math., Univ. of British Columbia, Vancouver, 1993. 23] J.C. Strikwerda. Finite Di erence Schemes and Partial Di erential Equations. Wadsworth & Brooks/Cole, 1989.

REFERENCES

27

24] L.N. Trefethen. Lax-stability vs. eigenvalue stability of spectral methods. In K.W. Morton and M.J. Bains, editors, Numerical Methods for Fluid Dynamics III, pages 237{253. Clarendon Press, 1988. 25] J.M. Varah. Stability restrictions on second order, three level nite di erence schemes for parabolic equations. Siam J. Numerical Analysis, 17(2):300{309, 1980. 26] R. Varga. Matrix Iterative Analysis. Prentice-Hall, Englewood Cli s, NJ, 1962.


赞助商链接
相关文章:
更多相关文章: