当前位置:首页 >> >>

Implicit-explicit Runge-Kutta methods for time-dependent partial differential equations


Implicit-Explicit Runge-Kutta Methods for Time-Dependent Partial Di erential Equations
Uri M. Ascher Steven J. Ruuthy Raymond J. Spiteriz March 13, 1997
Abstract

Implicit-explicit (IMEX) linear multistep time-discretization schemes for partial di erential equations have proved useful in many applications. However, they tend to have undesirable time-step restrictions when applied to convection-di usion problems, unless di usion strongly dominates and an appropriate BDF-based scheme is selected 2]. In this paper, we develop Runge-Kutta-based IMEX schemes that have better stability regions than the best known IMEX multistep schemes over a wide parameter range.

1 Introduction
When a time-dependent partial di erential equation (PDE) involves terms of di erent types, it is a natural idea to employ di erent discretizations for them. Implicit-explicit (IMEX) time-discretization schemes are an example of such a strategy. Linear multistep IMEX schemes have been used by many researchers, especially in conjunction with spectral methods 10, 3]. Some schemes of this type were proposed and analyzed as far back as the late 1970's 15, 5]. Instances of these methods have been successfully applied to the incompressible Navier-Stokes equations 9] and in environmental modeling studies 16]. A systematic, comparative study for PDEs of convection-di usion
Institute of Applied Mathematics and Department of Computer Science, University of British Columbia, Vancouver, BC, V6T 1Z4, Canada. (ascher@cs.ubc.ca). The work of this author was partially supported under NSERC Canada Grant OGP0004306. y Department of Mathematics, University of California, Los Angeles. (ruuth@math.ucla.edu). The work of this author was partially supported by an NSERC Postdoctoral Scholarship and NSF DMS94-04942. z Institute of Applied Mathematics and Department of Mathematics, University of British Columbia, Vancouver, BC, V6T 1Z2, Canada.(spiteri@math.ubc.ca). The work of this author was partially supported under NSERC Canada Grant OGP0004306.

1

1 INTRODUCTION

2

type was carried out in 2], and a corresponding study for reaction-di usion problems arising in morphology is reported in 11]. In this work, we consider problems of convection-di usion type (or hyperbolicparabolic equations), e.g.

ut = uux +

u;

> 0:

(1.1)

Discretization of the spatial derivatives (e.g. by nite-di erence, nite-element, nitevolume, or spectral methods) yields a very large ODE system in time,

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

(1.2)

where f corresponds to the convection (hyperbolic) term, uux, and g corresponds to the di usion (parabolic) term, u. An IMEX scheme consists of applying an implicit discretization for g and an explicit one for f . This is natural because the system (1.2) with f 0 is generally sti and linear, whereas with g 0 the system is not too sti and it is often nonlinear. Moreover, if one insists on using an implicit discretization for the latter case, the construction of iterative solvers is made challenging by the properties of the matrix to be inverted. In 2] we have analyzed and experimented with linear multistep schemes for (1.2). For di usion-dominated problems, schemes utilizing Backward Di erentiation Formulae (BDF) for g have optimal damping properties. The corresponding IMEX scheme is a semi-explicit BDF (SBDF) 15, 5, 9]. However, when di usion is not dominant, the choice of method is less obvious. Ideally, one would like a dissipative scheme for the hyperbolic term, so that the resulting IMEX scheme would have good stability and smoothing properties, independent of g(u). Yet, it is well known that this necessitates looking for schemes of order (or more accurately, backward steps) at least 3 (cf. 2]). The main di culty plaguing multistep IMEX schemes is the possibility that relatively small time steps become necessary due to stability restrictions (e.g. 14]). Such restrictions become worse when working with higher-order multistep schemes, as demonstrated in 2]. On the other hand, the stability region of explicit Runge-Kutta schemes actually increases slightly when three- and four-stage schemes are contemplated. Thus, we are led to investigate the development and performance of IMEX Runge-Kutta schemes. When the system (1.2) arises from a PDE in more than one spatial variable, the simplicity and e ciency of solving the algebraic equations corresponding to the implicit part of the discretization at each time step is of paramount importance. In 2], we show that, under appropriate circumstances and for certain IMEX schemes, not much more than one multigrid W-cycle per time-step is needed. Here, we aim at such a performance per Runge-Kutta stage. Hence, we concentrate on diagonallyimplicit Runge-Kutta (DIRK) schemes for g(u), paying particular attention to their attenuation properties. We construct a number of IMEX schemes of this sort, and investigate their properties. Since the e ective order of DIRK schemes drops to 1 for

2 IMEX RUNGE-KUTTA SCHEMES

3

very sti ODEs 8], we do not expect the Runge-Kutta schemes to be competitive with SBDF in the limit where g(u) is heavily dominant and very sti . However, our Runge-Kutta schemes are shown to have excellent properties for other parameter ranges. It is well known that, despite its recent rise in popularity, the use of Runge-Kutta schemes for the integration of time-dependent PDEs is not without drawbacks. One problem is a loss of accuracy order when time-dependent in ow boundary conditions are prescribed 4]. A remedy is proposed in 1]. Note also that, in the (usual) case where the accuracy of the solution at the end of the step is higher than at internal stages, the accuracy to which these internal stages must be calculated (say by an iterative method, as in x5) cannot be lowered. The system (1.2) can be cast as a partitioned system, v = f^(v; w) _ (1.3a) w = g(v; w) _ ^ (1.3b) with u = v + w, f^ = f (v + w) and g = g(v + w). The IMEX Runge-Kutta schemes ^ devised below are then partitioned Runge-Kutta schemes. Their orders are discussed in 7], xII.15. They can also be seen as a particular class of splitting methods. We note that for (1.2), v and w need not be known individually (indeed we may not have initial conditions for them). All we need to have is their smooth, local existence. In x2, we develop a number of IMEX Runge-Kutta schemes of up to four stages, which are up to third-order accurate. The schemes divide naturally into two classes { those whose last internal stage is identi ed with the solution at the next time instance (i.e. at the end of the time-step), and those where an additional quadrature is used at the end of the step. The rst class is particularly good for highly-sti problems, but it is the second class which seems to yield some of the more promising variants for a wide range of parameters, particularly the scheme identi ed as (3,4,3) (3 internal stages for the implicit formula, 4 stages for the explicit and a combined accuracy of order 3). In x3, we investigate the performance of the various schemes for a simple test equation, which arises from a von Neumann analysis of an advection-di usion equation. Stability regions are plotted in Figures 3.1 and 3.2. Results of test runs for a linear advection-di usion equation are reported in x4.1 and the Burgers equation is experimented with in x4.2. In x5, we investigate the performance of our various schemes in conjunction with employing a multigrid method for resolving the implicitness at each stage.

2 IMEX Runge-Kutta schemes
We now develop some IMEX Runge-Kutta schemes. For g, we consider an implicit s-stage DIRK scheme 8] with coe cients A 2 Rs s , c; b 2 Rs , in the usual Butcher

2 IMEX RUNGE-KUTTA SCHEMES

4

notation. Let = s + 1. For f , we consider an (s + 1)-stage explicit scheme with the ^ abscissae c = 0 and coe cients A 2 R ; ^ 2 R . To cast the DIRK scheme as ^ c b an (s +1)-stage scheme as well, we can pad the s-stage scheme with zeroes, obtaining the tableau 0 0 0 0 0 0 c 0 a 0 c 0 a a 0 . . . . ... . . . . . . . . . . . cs 0 as as ass 0 b b bs ~ Referring to the coe cients of this padded tableau as A 2 R ; ~ 2 R ; c 2 R , we b ~ see that c = c. This property simpli es signi cantly the form of the order conditions ~ ^ on the method's coe cients due to coupling between the individual schemes. One step from tn? to tn = tn? + k of the IMEX scheme is given as follows: Set ^ K = f (un? ): (2.1a)
1 11 21 2 22 1 2 1 2 1 1 1 1

For i = 1; : : : ; s do: Solve for Ki.

Ki = g(ui);
where

(2.1b)
i X j =1

ui = un? + k
1

i X j =1

ai;j Kj + k

a ^ ^i ;j Kj :
+1

(2.1c)

Evaluate ^ Ki Finally, evaluate
+1

= f (ui ):

(2.1d)

un = un? + k
1

s X j =1

bj Kj + k

X^ ^ bj Kj :
j =1

(2.1e)

We consider two special cases, which lead to two sub-families of methods: those which satisfy (2.2) and those which satisfy (2.3),(2.4):

2 IMEX RUNGE-KUTTA SCHEMES
1. In the case that ^ = ~ (and in particular ^ = 0), we have in place of (2.1e) b b b
1

5
s X j =1

un = un? + k
1 +1 +1

^ bj (Kj + Kj ):
+1

(2.2)

^ 2. In the case that ^s = 0, Ks need not be evaluated. Furthermore, if b bj = as;j ; ^j = as ;j ; j = 1; : : : ; s b ^ (2.3)
+1

(implying that the implicit scheme is sti y accurate and cs = 1), then substituting (2.1e) into the expression for us, we see that

un = us:
1 1

(2.4)

This is useful for very sti problems. We note also that the explicit scheme can now be generated from a general explicit s-stage scheme based on the abscissae c ; : : : ; cs? ,

c c

0
1 2

cs?

. . .

a ^ a ^
1

0 0 0 0 0 0 ^ a 0 0 . . . ... . . . . . . . . . as as as ^ ^ ^ 0 ^ ^ ^ ^s b b b b
21 31 32 1 2 3 1 2 3

0

In (2.1b), there are s nonlinear systems to solve involving g alone, and the operators to be inverted are Assuming that the Jacobian gu is evaluated only once (at tn? ), this operator is the same for each i if aii is independent of i, i.e. for SDIRK schemes. However, if an iterative method like 1-cycle multigrid is used for the implicit scheme, then perhaps it is less critical for the overall e ciency to insist on SDIRK, and a more general DIRK scheme can be contemplated. See 8], xIV.6, for sti y-accurate SDIRK schemes. In our search for particularly accurate and stable schemes, we have allowed the diagonal elements of A to di er from one another. Despite this, the schemes we recommend are all SDIRK, since there were su cient excess degrees of freedom in the solution process to allow for this computationally-e cient choice.
1

I ? kaiigu :

Example 2.1 Consider the case of (1.2) where u = q , f = f0 , g = f0 . p
1

Thus, we apply the implicit scheme to advance q and the explicit one to advance p. Observe that if f2 is in uenced by a control function, for instance, then at each stage

2

2 IMEX RUNGE-KUTTA SCHEMES
1

6

^ i, this in uence is propagated to Ki from Ki . This does not happen when using a \normal" explicit scheme (cf. 13]). Moreover, if f is independent of q then I ? kaiigu = I ?kaii(f )p 0 I
1

and the scheme becomes explicit! In particular, for Hamiltonian systems we can set u = (q; p)T , f = (0; ?Hq )T , g = (Hp; 0)T . If Hp is independent of q, then we obtain an explicit scheme which still retains some advantages of the implicit one.

2

We now proceed to construct some instances of our IMEX RK family of schemes. We will use the triplet (s; ; p) to identify a scheme (2.1), where s is the number of stages of the implicit scheme, is the number of explicit scheme stages (so = s + 1 for Case 1, i.e. (2.2), and = s for Case 2, i.e. (2.3)-(2.4) above), and p is the combined order of the scheme. The pair of backward and forward Euler schemes 1 1 0 0 1 1 can be padded to read 0 0 0 0 0 0 1 0 1 1 1 0 0 1 1 0 This yields the linear one-step IMEX (cf. 2]): un = un? + k(f (un? ) + g(un)): Note that this scheme satis es (2.3).
1 1

2.1 Forward-Backward Euler (1,1,1)

The above form of Euler's method does not satisfy ^ = ~. The variant that does is b b 0 0 0 0 0 0 1 0 1 1 1 0 0 1 0 1 i.e. un = un? + k(g(u ) + f (u )) where u = un? + kg(u ) + kf (un? ). This scheme involves an additional evaluation of f per step, compared to the previous one.
1 1 1 1 1 1 1

2.2 Forward-Backward Euler (1,2,1)

2 IMEX RUNGE-KUTTA SCHEMES

7

2.3 Implicit-Explicit Midpoint (1,2,2)
0 0 0 0 0 1
1 2 1 2

The Euler pairs are rst-order accurate and the (1,1,1) scheme has well-known disadvantages when g = 0 (cf. x3 and 2]). The following pair of implicit-explicit schemes 0 0 0 0 0 1
1 2 1 2

correspond to applying explicit midpoint for f and implicit midpoint for g. It is second-order accurate because the two schemes from which it is composed are each second-order accurate and c = c 7]. This scheme performs usually comparably to the ^ ~ unreasonably popular CNAB considered in 2], and has better symmetry properties.

Example 2.2 Applying the scheme (1,2,2) to the separable Hamiltonian system q_ = Hp(p); p = ?Hq (q) _
and de ning we obtain

q u = p ; f = Hp ; g = ?0 ; 0 Hq
^ K = Hp(pn? ) ; K = ?H0(q ) ; 0 q
1 1 1 1

q = qn? + k (K + K ) = qn? + k Hp(pn? ) ^ pn? pn? p 2 2 ?Hq (q )
1 1 1 1 1 1 1 1 1 1

q q ( qn? Hp ^ ^ K = Hp0p ) ; pn = pn? + k(K + K ) = pn? + k ?H(pq)) : n n? q(
2 1 1 1 1 2 1 1 1 1

We note that this scheme is explicit, and it can also be shown to be symplectic. It is identical to the leapfrog/Verlet scheme 12]

pn = pn? ? kHq (qn? + k=2Hp (pn? )) qn = qn? + kHp pn +2pn? :
1 1 1 1 1

2

2 IMEX RUNGE-KUTTA SCHEMES

8

The two-stage, third-order DIRK scheme with the best damping properties turns out to be the DIRK scheme (p. 207, 7]) 1?
3+ 3 6

2.4 A third-order combination (2,3,3)
1?2 1=2 0 1=2

with = p . For the test equation u = u, if we let z = k, then as z ! ?1, _ R(1) = 1 ? 3 ?0:7321. This is a marked improvement over the midpoint scheme (1,2,2) which has no attenuation at the sti ness limit 1. The corresponding third-order explicit Runge-Kutta scheme (ERK) is 0 1? 0 0 ? 1 2(1 ? ) 0 0 1=2 1=2 0 0 0

p

The resulting IMEX combination is third-order accurate, has some dissipation near the imaginary axis (like all three-stage third-order ERK methods) and has some attenuation of the stability function at 1. One may be concerned that the attenuation that the above scheme has may not be su cient for some problems. A two-stage, second-order DIRK scheme which is sti y accurate is (p. 106 8]) 1 1? 1? 0 0 1 0 0 0

2.5 L-stable, two-stage, second-order DIRK (2,3,2)

where =

2

?

p

2

2

. The corresponding three-stage second-order ERK is 1? 0 1? 0 0 0

Varying to get a dissipative stability region, we match the terms of the exponential in the stability function up to third order. This yields the stability region of a p three-stage, third-order explicit RK scheme, with = ?2 2=3. The resulting IMEX combination is second-order accurate.

2 IMEX RUNGE-KUTTA SCHEMES

9

We use the same form (with the same but with not yet speci ed) as in the previous IMEX scheme, except for requiring that (2.3) hold instead of (2.2), i.e., ^ = ; ^ = 1 ? ; ^ = 0. This gives a second-order scheme b b b 0 0 0 0 0 0 0 1 1? 1 1? 0 1? 1? 0 with = 1 ? .
1 2 3 1 2

2.6 L-stable, two-stage, second-order DIRK (2,2,2)

We would now like to exploit the larger dissipativity region, which four-stage, fourthorder ERK schemes have near the imaginary axis . A three-stage, third-order DIRK scheme which is sti y accurate is (p. 106 8]) 0 0 ? 0 1 b( ) b( ) b( ) b( ) where is the middle root of 6x ? 18x + 9x ? 1 = 0, b ( ) = ? + 4 ? and b ( ) = ? 5 + . Numerically, this evaluates to 0 0 :4358665215 :4358665215 :7179332608 :2820667392 :4358665215 0 1 1:208496649 ?:644363171 :4358665215 1:208496649 ?:644363171 :4358665215 The corresponding four-stage, third-order ERK, constructed such that it has the same (desirable) stability region as all four-stage, fourth-order ERK schemes, is 0 0 0 0 0 0 0 0 ^ a () a () 0 0 ^ 1 a () a ^ ^ a 0 ^ 0 b( ) b( ) The third-order conditions indicate that this is a two-parameter family of schemes. Taking our degrees of freedom to be a ; a , the remaining expressions are ^ ^
(1+ ) 2 (1 1 ) 2 2 1 2 3 2 2 3 2 2 5 4 1 3 2 2 1 4 (1+ ) 2 31 32 41 42 43 1 2 42 43

2.7 L-stable, three-stage, third-order DIRK (3,4,3)

3 a = 1? 9 + 2 ^ 2
31

2

a + 11 ? 21 + 15 ^ 4 2 4
42

2

9 a ? 7 + 13 ? 2 ; ^ 2
43 2

2 IMEX RUNGE-KUTTA SCHEMES a = ?1 + 9 ? 3 ^ 2 2
32 2

10

a ^ + ? 11 + 21 ? 15 4 2 4
42 41 42 43

2

a + 4 ? 25 + 9 ; ^ 2 2
43 2

exponential in the stability function up to fourth order, we obtain the scheme 0 0 0 0 0 :4358665215 :4358665215 0 0 0 :7179332608 :3212788860 :3966543747 0 0 1 ?:105858296 :5529291479 :5529291479 0 0 1:208496649 ?:644363171 :4358665215 The resulting IMEX combination is third-order accurate.

a =1?^ ?^ : ^ a a Selecting particular values for ^ ( ) and ^ ( ) so as to match the terms of the a a
42 43

2.8 A four-stage, third-order combination (4,4,3)

No pair consisting of a three-stage, L-stable DIRK and a four-stage ERK satisfying (2.3) with a combined third-order accuracy could be found. In order to obtain a third-order formula satisfying (2.3), we need to go to a four-stage, third-order, Lstable DIRK, coupled with a four-stage ERK in the manner indicated in Case 2 of x2. After having satis ed the order conditions and imposed L-stability, we are left with a few families of schemes. Experimenting with the test equation of x3, we have selected a scheme which has rational coe cients which are not too large, and a diagonal entry which is close to that of the four-stage, fourth-order, L-stable DIRK. In particular, we have chosen ^ c = , accompanied by the choices c = , b = , and a = . The resulting scheme is
1 1 2 3 1 2 4 1 2 43 1 2 1 2 2 3 1 2

1 0 0 1
1 2 2 3 1 2

?

1 2 1 6 3 2 3 2

1 2

? ?

0 0 0 0 0 0
1 2 1 2 3 2 3 2 1 2 1 2 1 2 1 2 1 2

1 2 11 18 5 6 1 4 1 4

0 0 0 0 0 0 0 0 ? 0
1 18 5 6 7 4 7 4 1 2 3 4 3 4

0 0 0 0 ? 0 ? 0
7 4 7 4

3 TEST EQUATION FOR CONVECTION-DIFFUSION

11

3 Test equation for convection-di usion
As is shown in 2], using a centered discretization scheme for the spatial derivatives of advection-di usion PDEs, a von Neumann analysis yields that a simple yet appropriate test equation is the scalar ODE (1.2) with

f = i u; g = u; p where ; are real constants, usually 0, > 0, and i = ?1. For a given time step k, we can de ne x = k ; y = k ; z = x + iy;
and write an IMEX step (2.1)-(2.2) as

1

un = un? + z
1

s X j =1

bj uj R(z)un? ;
1

(3.1)

where (1 + iy^i a
+1;1 1

)un? + ij? (xaij + iy^i ;j )uj a i = 1; : : : ; s: (3.2) ui = 1 ? xaii For schemes where = s, i.e. where (2.3) holds, an expression like (3.1) need not be used because un = us. Let us consider the two Euler schemes rst. For the usual pair (1,1,1), equation (3.2) implies y un = 1 + ix un? : 1?
1 =1 +1 +1 1 2 1

P

On the imaginary axis, junj (1 + y )jun? j, and the scheme is unconditionally unstable for any y 6= 0 (cf. 2]). As x ! ?1 and with y bounded, however, the scheme is unconditionally stable and un ! 0, i.e. attenuation compatible with that of the exact solution is realized. For the other pair, (1,2,1), we get from (3.1)-(3.2)

y un = 1 + z 1 + ix un? : 1?
1

Setting x = 0 initially, we obtain that, along the imaginary axis,

junj
1

2

(1 ? y + y )jun? j ;
2 4 1 2

W. Hunsdorfer (Workshop on Innovative Time Integrators, CWI, Amsterdam, 1996) has indicated that a more general test equation where and are allowed to be complex is useful, but we stay with the simpler test equations and lend further support to the analysis with more general numerical experiments in the following sections.

3 TEST EQUATION FOR CONVECTION-DIFFUSION
10 8
Time-Step Restrictions

12

k

6 4 2

(1,1,1) (2,2,2) (4,4,3)

?10

0

2

?10

1

?10
=

0

?10?

1

?10?

2

Figure 3.1: Time-step stability restrictions for various values of and , for IMEX schemes satisfying (2.3). Note that all these schemes perform well in the highattenuation limit. so stability is achieved provided that

jyj 1:
This corresponds to the CFL condition. On the other hand, letting x ! ?1, we also obtain the restriction jyj 1 in the high-attenuation limit. Thus, we see that the variant (1,2,1) has stability characteristics which are preferable over the usual backward-forward Euler scheme in the hyperbolic limit and worse in the parabolic limit. Now, we can ask, given a xed ratio x=y = = , what is the maximum y > 0 for which the scheme is stable, i.e. jR(z)j < 1? This gives the stability restriction on the step-size for a given problem. In Figure 3.1, we plot the resulting curves for the IMEX scheme presented earlier of the second type, i.e. those satisfying (2.3). Corresponding curves for the schemes satisfying (2.2) are plotted in Figure 3.2. Common to all the schemes in Figure 3.2 is that, as ! ?1, there is a time-step restriction involving , typically k < 1. In contrast, no such restriction occurs for schemes satisfying (2.3), and hence (2.4), in Figure 3.1. The behavior of these latter schemes is qualitatively similar in the high-attenuation limit to that of the SBDF schemes 2]. On the other hand, the scheme (3,4,3) in Figure 3.2 suggests a potential advantage, in terms of allowable time steps, over a large range of parameter values.

4 FINITE-DIFFERENCE APPROXIMATIONS IN 1D
12 10 8
Time-Step Restrictions (1,2,1) (1,2,2) (2,3,3) (2,3,2) (3,4,3)

13

k

6 4 2 0

?10

2

?10

1

?10
=

0

?10?

1

?10?

2

Figure 3.2: Time-step stability restrictions for various values of and , for IMEX schemes satisfying (2.2). Note that some of these schemes, especially (3,4,3), perform well when ? is not too large in magnitude, but they all have time-step restrictions in the high-attenuation limit.

4 Finite-Di erence Approximations in 1D
Consider the one-dimensional variable-coe cient problem ut + sin(2 x)ux = uxx; (4.1) subject to periodic boundary conditions on the interval 0,1] and initial condition u(x; 0) = sin(2 x): We use centered, second-order di erences for the spatial derivatives. Note that the solution is smooth for all 0.

4.1 A linear problem

Large-viscosity examples
2

To examine the behavior of IMEX RK schemes for large viscosities , the model problem (4.1) was approximated using the time step k = 1:8h , h = . For the spatial
2 0 0 1 63

By large viscosity, we mean that the mesh Reynolds number is small (< 2; see e.g. 2]). An estimate for the mesh Reynolds number is given by R = ah where represents viscosity, a represents characteristic speed, and h represents grid spacing.

4 FINITE-DIFFERENCE APPROXIMATIONS IN 1D
0 0 0 0

14

discretization, we chose a grid spacing of h = h ; h =2; h =4 or h =8. Utilizing several schemes, computations to time t = 2 were performed for viscosities, , in the range :01 :1. For h = h , these values correspond to mesh Reynolds numbers, R, in the range 1:59 R 0:159. For smaller h, R is correspondingly smaller. The relative errors measured in the maximum norm are plotted against in Figure 4.1. In Figure 4.1(a), the errors for the three- and four-stage schemes are omitted since they are always within a factor of 2 of the two-stage, second-order scheme (2,3,2). As expected from the theory of the previous section, the three-stage DIRK (3,4,3) and the schemes satisfying property (2.3) allow for the largest stable time steps when is large. When we decrease h, keeping k xed, the methods (1,2,2), (2,3,3) and (2,3,2) become unstable over the entire depicted range of . Among the stable methods, the higher order of (3,4,3) and (4,4,3) becomes more apparent as the spatial error ceases to dominate the temporal error. For h = h =8, only the methods satisfying (2.3) remain stable. No e ective error reduction is detected for these values of h and k. Doubling the time step while keeping h = h also causes the two-stage methods (2,3,2) and (2,3,3) to become unstable for almost the entire -interval depicted. Nonetheless, the three- and four-stage methods remain stable over the entire interval and the scheme (2,2,2) is stable over most of the interval. Even upon tripling the time step (see Figure 4.2), the three- and four-stage schemes are appropriate for approximating strongly-damped ows. Indeed, a comparison with the results of 2] indicates that the three-stage method (3,4,3) is more accurate and stable in this example than the popular CNAB scheme with one-third the step size, even for this relatively coarse step size k.
0 0 0

To examine the behavior of IMEX RK methods for large mesh Reynolds numbers, example (4.1) was approximated using discretization step sizes h = in space and k = 1:8h 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. The relative errors in maximum norm for these computations are plotted against in Figure 4.3. Here, we nd that the errors for the third-order schemes and the two-stage, second-order DIRK (2,3,2) nearly coincide since they arise predominantly from the spatial discretization. It is especially noteworthy that the errors for these four schemes are always smaller than those arising for the two-step IMEX schemes considered in 2] even though twice the time step is applied.
1 81

Small-viscosity examples

5 TIME-DEPENDENT MULTIGRID IN 2D

15

4.2 Burgers Equation

Qualitatively similar results to those obtained above for the linear problem (4.1) are obtained for the Burgers equation

ut + uux = uxx; u(x; 0) = sin( x);
0 < x < 1:

(4.2)

with which we have experimented, subject to homogeneous Dirichlet boundary conditions and the initial conditions For small , a boundary layer develops near x = 1. This necessitates the use of upwinding. Note that upwinding can be viewed as an addition of a discrete O(h) di usion term (where h is the spatial step-size) to a centered-di erence approximation of uux.

5 Time-dependent Multigrid in 2D
We consider the two-dimensional convection-di usion problem,

u + (u r)u =
t

u;

(5.1) 0,1] 0,1] and

where u (u; v). We carry out our computations on the square consider periodic boundary conditions and the initial conditions

u(x; y; 0) = sin 2 (x + y)] + 0:005 cos 2 (64x + 63y)] v(x; y; 0) = sin 2 (x + y)] + 0:005 cos 2 (64x + 63y)]:
For the spatial discretization, we use standard second-order centered di erences. Time-stepping is carried out using a variety of IMEX RK 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 be solved at each stage. Such systems are solved e ciently using a multigrid algorithm, the components of which are outlined in 2]. The model problem (5.1) was approximated using a mesh size h = and a residual tolerance of TOL=0.003 at each stage. The time step was selected to be k = 0:00625 to achieve stable results. For several IMEX RK schemes, the average number of ne-grid iterations at each time step was computed. The results for V(1,1)cycles are provided in the table below. In addition to the schemes developed in x2, we include results for Griepentrog's scheme 6, 7], which is not of the form (2.1).
1 128

REFERENCES
Scheme Forward-backward Euler (1,1,1) Forward-backward Euler (1,2,1) Implicit-explicit midpoint (1,2,2) Third-order combination (2,3,3) Griepentrog's scheme L-stable, 2-stage, 2nd-order DIRK (2,3,2) L-stable, 2-stage, 2nd-order DIRK (2,2,2) L-stable, 3-stage, DIRK (3,4,3) L-stable, 4-stage, 3rd-order scheme (4,4,3)

16 0.01 1.94 1.97 2.23 4.32 4.58 3.81 3.90 5.87 7.84
viscosity, 0.02 0.03 0.04 1.90 1.71 1.29 1.94 1.48 1.19 2.87 3.13 3.29 4.42 3.68 3.71 5.23 5.39 5.61 3.35 2.74 2.58 3.35 2.81 2.58 5.35 3.81 3.61 7.35 5.16 4.81

0.05 1.19 1.19 3.68 3.81 5.90 2.61 2.58 3.58 4.77

From the table, we see that the strongly-damping L-stable schemes require the fewest ne-grid iterations per stage for large-viscosity problems. Weakly-damping schemes require far more e ort to solve the implicit equations accurately, because lingering high-frequency modes necessitate more work on the nest grid. This is especially evident for the implicit-explicit midpoint scheme (1,2,2) since more than three iterations were required per stage when 0:03. Griepentrog's scheme also requires extra ne-grid iterations, presumably to damp out the high-frequency components which arise during the initial stage.

References
1] S. Abarbanel, D. Gottlieb, and M. Carpenter. On the removal of boundary errors caused by Runge-Kutta integration of nonlinear partial di erential equations. SIAM J. Scient. Comput., 17:777{782, 1996. 2] U. Ascher, S. Ruuth, and B. Wetton. Implicit-explicit methods for timedependent PDE's. SIAM J. Numer. Anal., 32:797{823, 1995. 3] C. Canuto, M.Y. Hussaini, A. Quarteroni, and T.A. Zang. Spectral Methods in Fluid Dynamics. Springer-Verlag, 1987. 4] M. Carpenter, D. Gottlieb, S. Abarbanel, and W. Don. The theoretical accuracy of Runge-Kutta time discretizations for the initial-boundary value problem: A study of the boundary error. SIAM J. Scient. Comput., 16:1241{1252, 1995. 5] M. Crouzeix. Une methode multipas implicite-explicite pour l'approximation des equations d'evolution paraboliques. Numer. Math., 35:257{276, 1980. 6] E. Griepentrog. Gemischte Runge-Kutta-verfahren fur steife systeme. In Seminarbereicht Sekt. Math., number 11, pages 19{29, Humboldt-Universitat Berlin, 1978.

REFERENCES

17

7] E. Hairer, S.P. Norsett, and G. Wanner. Solving Ordinary Di erential Equations I: Nonsti Problems. Springer-Verlag, 1993. 8] E. Hairer and G. Wanner. Solving Ordinary Di erential Equations II: Sti and Di erential-Algebraic Problems. Springer-Verlag, 1991. 9] 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. 10] J. Kim and P. Moin. Application of a fractional-step method to incompressible Navier-Stokes equations. J. Computational Phys., 59:308{323, 1985. 11] S. Ruuth. Implicit-explicit methods for reaction-di usion problems in pattern formation. J. Math. Biology, 34(2):148{176, 1995. 12] Robert D. Skeel, Guihua Zhang, and Tamar Schlick. A family of symplectic integrators: stability, accuracy, and molecular dynamics applications. SIAM J. Sci. Comput., 18(1):203{222, 1997. 13] R. Spiteri, U. Ascher, and D. Pai. Numerical solution of di erential systems with algebraic inequalities arising in robot programming. In Proceedings of the IEEE Conference on Robotics and Automation, 1995. 14] S. Turek. A comparative study of time stepping techniques for the incompressible Navier-Stokes equations: from fully implicit non-linear schemes to semi-explicit projection methods. Int. J. Num. Meth in Fluids, 22:987{1011, 1996. 15] 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. 16] J.G. Verwer, J.G. Blom, and W. Hundsdorfer. An implicit-explicit approach for atmospheric transport-chemistry problems. Applied Numerical Mathematics, 20:191{209, 1996.

REFERENCES

18

10

0

10 Implicit?explicit midpoint (1,2,2) Third?order combination (2,3,3)

0

10 Relative Error

?1

10 Relative Error

?1

(2,2,2)

2?stage, 2nd?order DIRK (2,2,2)

10

?2

10

?2

(3,4,3) (4,4,3)

2?stage, 2nd?order DIRK (2,3,2)

10

?3

10

?2

10 Viscosity

?1

10

?3

10

?2

10 Viscosity

?1

1 (a) h = 63

1 (b) h = 126

10

0

10

0

(3,4,3) 10
?1

10 (2,2,2)

?1

(2,2,2) Relative Error Relative Error

10

?2

10

?2

(4,4,3)

(4,4,3)
?3 ?3

10

10 (3,4,3)

10

?4

10

?2

10 Viscosity

?1

10

?4

10

?2

10 Viscosity

?1

1 (c) h = 252

1 (d) h = 504
1 63

Figure 4.1: Large-viscosity behavior for k = 1:8h ; h =
0 0

REFERENCES
10
0

19

2?stage, 2nd?order DIRK (2,2,2)

10 Relative Error

?1

4?stage, 3rd?order scheme (4,4,3)

10

?2

3?stage DIRK (3,4,3)

10

?3

10

?2

10 Viscosity

?1

Figure 4.2: Large-viscosity behavior for k = 5:4h; h = h
10
0

0

Forward?backward Euler (1,1,1), (1,2,1)

10 Relative Error

?1

Implicit?explicit midpoint (1,2,2) and 2?stage, 2nd?order DIRK (2,2,2)

Multistage schemes: (2,3,3), (2,3,2), (3,4,3), (4,4,3)

10

?2

10

?3

10

?3

10 Viscosity

?2

Figure 4.3: Small-viscosity example for k = 1:8h; h =

1 81


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