IMPLICIT-EXPLICIT METHODS FOR TIME-DEPENDENT PDE'SBySteven J. RuuthBMATH, University of Waterloo, 1991A THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFMASTER OF SCIENCEinTHE FACULTY OF GRADUATE STUDIESDEPARTMENT OF MATHEMATICSANDINSTITUTE OF APPLIED MATHEMATICSWe accept this thesis as conformingto the required standardTH UN ERSITY OF BRITISH COLUMBIAApril 1993© Steven J. Ruuth, 1993In presenting this thesis in partial fulfilment of the requirements for an advanced degree atthe University of British Columbia, I agree that the Library shall make it freely availablefor reference and study. I further agree that permission for extensive copying of thisthesis for scholarly purposes may be granted by the head of my department or by hisor her representatives. It is understood that copying or publication of this thesis forfinancial gain shall not be allowed without my written permission.Department of MathematicsThe University of British Columbia2075 Wesbrook PlaceVancouver, CanadaV6T 1W5Date: AbstractVarious methods have been proposed to integrate dynamical systems arising from spa-tially discretized time-dependent partial differential equations. For problems with termsof different types, implicit-explicit (IMEX) schemes have been used, especially in con-junction with spectral methods. For convection-diffusion problems, for example, onewould use an explicit scheme for the convection term and an implicit scheme for thediffusion term. Reaction-diffusion problems can also be approximated in this manner. Inthis work we systematically analyze the performance of such schemes, propose improvednew schemes and pay particular attention to their relative performance in the context offast multigrid algorithms and aliasing reduction for spectral methods.For the prototype linear advection-diffusion equation, a stability analysis for first,second, third and fourth order multistep IMEX schemes is performed. Stable schemespermitting large time steps for a wide variety of problems and yielding appropriate decayof high frequency error modes are identified.Numerical experiments demonstrate that weak decay of high frequency modes canlead to extra iterations on the finest grid when using multigrid computations with finitedifference spatial discretization, and to aliasing when using spectral collocation for spa-tial discretization. When this behaviour occurs, use of weakly damping schemes suchas the popular combination of Crank-Nicolson with second order Adams-Bashforth isdiscouraged and better alternatives are proposed.Our findings are demonstrated on several examples.iiTable of ContentsAbstract^ iiList of Figures viAcknowledgements^ viii1 Introduction 11.1 IMEX Schemes ^ 11.2 Multigrid and IMEX Schemes ^ 31.3 Spectral Methods and IMEX Schemes 51.4 Overview^ 62 First and Second Order IMEX Schemes 82.1 General Linear Multistep IMEX Schemes ^ 82.2 Optimal Accuracy for IMEX Schemes 92.3 Method of Analysis ^ 122.4 Schur and von Neumann Polynomials ^ 142.5 First Order IMEX Schemes ^ 152.6 Stability of First Order IMEX Schemes ^ 162.7 Second Order IMEX Schemes ^ 172.8 Choice of c for Second Order Schemes 192.9 Second Order Instability along the Nonzero i3-Axis ^ 212.10 Second Order Stability Contours ^ 22iii32.11 Choice of Second Order Method ^Higher Order IMEX Schemes24293.1 Third Order IMEX Schemes ^ 293.2 Choice of parameters 303.3 Third Order Stability along the 1 3-Axis ^ 323.4 Third Order Stability Contours 333.5 Fourth Order IMEX Schemes ^ 383.6 Comparison to Lower Order Methods 404 Convection-Diffusion Experiments 424.1 Finite Difference Approximations ^ 424.1.1^Variable Coefficient Problem 424.1.2^Small Mesh Reynolds Number Computations ^ 434.1.3^Large Mesh Reynolds Number Computations 444.1.4^Zero Viscosity Computations ^ 464.2 Spectral Collocation ^ 484.2.1^The Burgers Equation ^ 484.2.2^Fourier Spectral Collocation 514.2.3^Chebyshev Spectral Collocation ^ 534.2.4^Aliasing in Pseudospectral Methods 534.3 Time-Dependent Multigrid ^ 574.3.1^The Convection-Diffusion Problem ^ 574.3.2^Components of the Multigrid Algorithm 584.3.3^Relative Efficiencies of Several IMEX Schemes ^ 594.4 Choice of IMEX Scheme ^ 604.4.1^Finite Differences 60iv4.4.2 Spectral Methods ^625 Reaction-Diffusion Experiments^ 645.1 Reaction-Diffusion Systems ^645.2 Treatment of Reaction Terms ^ 665.3 Computations ^ 695.3.1 The Schnackenberg System ^ 695.3.2 One-Dimensional Computations 705.3.3 Two-Dimensional Computations ^ 745.4 Choice of IMEX Scheme ^ 766 Conclusions^ 796.1 Results ^796.2 Future Work^ 82Bibliography^ 84vList of Figures2.1 Half Ellipse of Possible (a, /9) ^132.2 Amplification over a Time Interval of Length k for x = i f3x + ax ^ 142.3 Stability Contours for y = 1 (First Order) ^172.4 Stability Contours for y = 0.5 (First Order) ^182.5 Stability Contours for CNLF, (7,c) = (0,1) ^222.6 Stability Contours for CNAB, (y, c) = (1., 0) ^232.7 Stability Contours for Modified CNAB, (7,c) = (2, ) ^232.8 Stability Contours for SBDF, (y, c) = (1,0) ^242.9 Stability Contours for (y, c) = (.25, .28125) ^252.10 Stability Contours for (y, c) = (.25, .75) ^252.11 Time Step Restriction For Various y where c = 1 — 7 ^ 272.12 Time Step Restriction For Various 7 where c = (1 7)2 ^ 282.13 Time Step Restriction For Various y where c = 0 283.14 Stability Contours for (y, 0, c) = olim (y, 0, .46610) ^353.15 Stability Contours for (-y,0,c) = (0, —2.036, —.876) ^353.16 Stability Contours for (y, 0, c) = (2, 0,0) ^363.17 Stability Contours for (y, 0, c) = (-, —1.21, —.5) ^363.18 Stability Contours for Third Order SBDF, (7,0,c) = (1,0,0) ^373.19 Zoom-in Around a = 0 for Third Order SBDF ^ 373.20 Stability Contours for Fourth Order SBDF 393.21 Zoom-in Around a = 0 for Fourth Order SBDF ^ 39vi3.22 Time Step Restrictions for Various IMEX Schemes ^ 414.23 Large Viscosity Behaviour For Various IMEX Methods ^ 454.24 Second Order Results for 7 = J -4 for Various c ^ 454.25 Small Viscosity Example For Various IMEX Methods ^ 474.26 Zero Viscosity Example For Various Second Order Methods, y ^ 494.27 Zero Viscosity Example For Higher Order Methods ^ 494.28 Burgers Equation for Various t for v = 0.01 504.29 Fourier Spectral Collocation for Burgers Equation ^ 524.30 Chebyshev Spectral Collocation for Burgers Equation 544.31 Aliased and Non-Aliased Computations for CNAB ^ 564.32 Multigrid V-Cycle Iterations, v = 0.03, h = -1218 ^604.33 Multigrid W-Cycle Iterations, v = 0.03, h = 1218 ^614.34 Multigrid W-Cycle Iterations, v = 0.02, h = 32 ^ 615.35 Stability Contours for Second Order Adams-Bashforth ^ 675.36 Stability Contours for the Explicit Extrapolation of SBDF ^ 685.37 Contours for the Explicit Extrapolation of 3rd Order SBDF ^ 685.38 Concentration, u, at t = 0.8 for v = 60 ^715.39 Concentration, v, at t = 0.8 for v = 60 ^715.40 Concentration, u, at t = 0.8 for v = 160 ^725.41 Concentration, v, at t = 0.8 for v = 160 ^725.42 Relative Errors for Several IMEX Schemes for Various v and h = 210^745.43 IMEX Computations for v, at t = 0.8 for v = 160 755.44 Relative Errors for Several IMEX Schemes for Various v and h = so 755.45 Concentration, u, at t = 0.2 for v = 120 77viiAcknowledgementsFirst and foremost, many thanks to my supervisors Dr. Uri Ascher and Dr. Brian Wettonfor their many helpful suggestions while working on this thesis. I would also like to thankNSERC for financially supporting me during my graduate work.viiiChapter 1Introduction1.1 IMEX SchemesThis thesis considers multistep methods for computing approximate solutions to largesystems of ordinary differential equations which arise from the spatial discretization oftime-dependent partial differential equations. Such systems typically have the formzi = f (u) + vg(u)^ (1.1)where II glj is normalized and v is a nonnegative parameter.The term f(u) in (1.1) is some possibly nonlinear term which we do not want tointegrate implicitly. This could be because the Jacobian of f(u) is non-symmetric andnon-definite and an iterative solution of the implicit equations is desired, or the Jacobiancould be dense, as in spectral methods, requiring the inversion of a full matrix at eachtime step. One may simply wish to integrate f (u) explicitly for ease of implementation.The term vg(u), however, is a stiff term which should be integrated implicitly to avoidexcessively small time steps. Frequently vg(u) is a linear diffusion term, in which casethe implicit equations form a linear system which is positive definite, symmetric andsparse. Such systems can be solved efficiently by iterative techniques (see Varga[23] orHackbusch[10]). Thus for problems of the form (1.1) it often makes sense to integratevg(u) implicitly and f (u) explicitly. Such implicit-explicit methods will be referred to asIMEX schemes.1Chapter 1 . Introduction^ 2Problems of the form (1.1) frequently occur in practice. For example, convection-diffusion problems lead to such systems. Kim and Moin[13], for instance, applied aprojection scheme to the incompressible Navier-Stokes equations to obtain a nonlinearconvection-diffusion problem to solve at each time step. This system was solved us-ing second order Adams-Bashforth for the convection term and Crank-Nicolson for thediffusion term. Applied to (1.1) this popular scheme givesn+1 3 1u — un = .-f (un) — .- f (un- ') + l [g (un+1 ) + g(un )]where k is the constant discretization step size and u" is the numerical approximation tou(kn).A wide variety of other applications for IMEX schemes are also possible. For example,solutions to reaction-diffusion systems arising in chemistry and mathematical biology canbe computed using this technique. For these problems the nonlinear reaction term canbe treated explicitly while the diffusion term is treated implicitly. Examples of reaction-diffusion systems from a biological standpoint can be found in Murray[16].Several authors have analyzed specific IMEX schemes. For example, an experimentalanalysis of several implicit-explicit schemes including (1.2) was carried out by Basdevantet al.[1]. Varah[22] determines some stability properties for certain second order IMEXschemes. These studies do not address how to choose the best IMEX scheme for agiven system (1.1). This thesis seeks efficient IMEX schemes for convection-diffusion andreaction-diffusion problems by systematically examining the stability properties of thelinearized systems and performing numerical experiments.Methods which allow the largest stable time steps are sought. In particular, forv > 1 we seek methods possessing a mild time step restriction since the system (1.1) isdominated by the implicitly handled "diffusion term". For other values of v > 0 we seekschemes which have reasonable time-stepping restrictions. A restriction similar to thek(1.2)Chapter 1. Introduction^ 3Courant-Friedrich-Lewy (CFL) condition is reasonable because this guarantees that thedomain of dependence of the scheme includes the domain of dependence of the differentialequation for explicit schemes applied to the one-dimensional wave equation, u t = aus(see Strikwerda[19]).Strong decay of high frequency spatial modes is another property which is importantin many applications. Schemes which possess this property are identified and comparedto standard schemes such as (1.2). As discussed in the next two sections, this property hasimportant implications for the efficiency of time-dependent multigrid and pseudospectralcomputations.1.2 Multigrid and IMEX SchemesAs mentioned in the previous section, IMEX schemes typically generate at each timestep implicit systems represented by sparse, symmetric, positive definite matrices. Intwo or higher dimensions, iterative techniques are particularly suitable for solving suchsystems. A multigrid solution for the implicit equations may be especially attractivebecause only 0(N) floating point operations are ideally needed for solving simple steadystate problems, where N is the size of the ODE system (1.1) arising from N spatial gridpoints (see Hackbusch[10]).In many physical problems with diffusion, high frequency solution components dissi-pate quickly. For multigrid this is a desirable property for the numerical solution since itmay reduce the number of costly iterations on the finest grid. Brandt and Greenwald[5]argue that non-physical high frequency components are especially harmful when applyingred-black Gauss-Seidel relaxation. This is because this smoothing technique aliases highfrequency errors into low frequency components which linger on for many time steps.Nonetheless, red-black Gauss-Seidel relaxation is often used, since it has a very goodChapter 1. Introduction^ 4smoothing rate and it is parallelizable.Several useful suggestions have been made by Brandt and Greenwald[5] to addressthese problems. One of these is use of a double discretization technique. This method usesBackward Euler during relaxation sweeps and Crank-Nicolson for calculating transferredresiduals. We do not consider this technique because it may be fragile in applicationssince high frequency components must be damped very strongly to obtain second orderaccuracy using a first order scheme for relaxation sweeps.Brandt and Greenwald[5] recommend use of W-cycles instead of V-cycles to help re-move aliased low frequency components from coarser grids. W-cycles are more expensivethan V-cycles, however. Furthermore, this method does nothing to remove non-physicalhigh frequencies on the finest grid.They further propose to avoid the final smoothing pass at each time step. This mayreduce low frequency aliasing errors since it ensures that all relaxations are followed bya coarse grid correction. However, this idea also does not eliminate non-physical highfrequency components on the finest grid.Finally, Brandt and Greenwald[5] suggest a modified FMG algorithm where the firstcoarse grid correction is performed before iterating on the finest grid. This should pro-duce the maximum effect from relaxation since the increment between time steps isgenerally smooth. In Chapter 4, however, numerical experiments for the 2D convection-diffusion problem demonstrate that this method alone can be inefficient for eliminatingnon-physical high frequency components.This thesis proposes use of a highly damping time-stepping method to control non-physical high frequency spatial components. This technique produces little or no addi-tional expense because strongly damping time-stepping schemes require approximatelythe same number of floating point operations as standard schemes such as (1.2). In Chap-ter 4, highly damping IMEX schemes are compared experimentally to standard schemesChapter 1. Introduction^ 5to verify that these methods are very effective for reducing the number of iterationson the finest grid, thus providing a more efficient multigrid method for time-dependentproblems.1.3 Spectral Methods and IMEX SchemesAs mentioned in Section 1.1, IMEX schemes are frequently used in spectral methods.Fully implicit schemes are impractical because a dense system originating from the con-vection term must be solved at each time step. Explicit schemes are usually avoidedbecause they frequently require unreasonably small time steps to ensure stability. Thetime step restriction for a Chebyshev basis (see Canuto et al.[6]) is particularly severe.For this basis, the one-dimensional heat equation with Dirichlet or Neumann boundaryconditions, for example, requires the time step restriction k = 0(N - 4 ) as N ---+ oo whereN is the number of basis functions. This result follows from the fact that the largesteigenvalue, A, has magnitude 1A1 = 0(N4 ) as N ---> oo (see Canuto et al.[6]). For furtherinformation about spectral methods applied to fluids and other problems see Canuto etal.[6] or Boyd[2].It is known that aliasing effects can complicate solutions for pseudospectral methods(see Canuto et al.[6]). These effects occur when nonlinear terms produce high frequenciesthat cannot be represented in the basis and thus contribute erroneously to low frequencies.Since only the highest frequencies alias into low frequencies we expect that a weak decayof high frequency spatial modes will contribute excessively to aliasing effects.Several methods have been proposed to deal with this behaviour. One such method isthe 3/2's rule (see Canuto et al.[6]). Using this technique, a quadratically nonlinear termis computed using 3N/2 basis functions. At the end of the calculation, coefficients for thebasis functions representing the N lowest frequencies are kept. Only 3N/2 basis functionsChapter 1. Introduction^ 6are required since frequencies which cannot be represented in the basis contribute tofrequencies that are later neglected. An alternative to the 3/2's rule is the 2/3's rule. Forthis method, the highest 1/3 frequencies are eliminated before any aliasing calculation.These methods eliminate aliasing but the 3/2's rule requires approximately 50% morefloating point operations than an aliased computation (see Canuto et al.[6]), while the2/3's rule causes a dramatic loss of high frequency information at each time step. Othermethods (see Canuto et al.[6]) such as applying a finer spatial mesh or using phase shiftsalso necessitate extra work.We propose to reduce aliasing by using a highly damping time-stepping method. InChapter 4, numerical experiments are performed to compare highly damping schemes tostandard schemes such as (1.2). These computations demonstrate that this inexpensivemethod of reducing aliasing can be very effective in aliased pseudospectral computations.1.4 OverviewIn Chapter 2, general linear multistep IMEX schemes are defined, and the optimal ac-curacy for such schemes is derived. This is followed by a description of characteristicpolynomials and other aspects of the stability theory used for analyzing properties ofmultistep schemes. First and second order schemes are discussed. In particular, parame-terizations for first and second order schemes are given. For the linear advection-diffusionequation, stability properties for several popular methods are determined and stabilityfor v = 0 is considered. High frequency modal decay and time step restrictions for sec-ond order schemes are studied in detail to determine the optimal method for variousproblems.In Chapter 3, third and fourth order IMEX schemes are discussed. A parameterizationfor third order schemes is given and the method with the strongest high frequency decayChapter 1. Introduction^ 7is identified. The stability properties for several methods are compared for v > 1 andv = 0. A fourth order scheme is defined and several of its stability properties are derivedfor v >> 1 and v = 0. This chapter concludes with a comparison of second, third andfourth order IMEX schemes.In Chapter 4, numerical experiments for convection-diffusion problems are carriedout. For the finite difference case, solutions to variable coefficient and nonlinear problemsare computed to support the stability results of Chapter 2. Some properties of IMEXschemes in pseudospectral methods are also demonstrated for the Burgers equation. Inparticular, these computations show that the use of strongly damping time-steppingschemes reduces aliasing inexpensively and effectively. By computing the solutions tothe 2D convection-diffusion equation using finite differences, it is also shown that the useof a strongly damping time-stepping scheme can give a more efficient multigrid methodfor time-dependent problems. This chapter concludes by discussing how to choose IMEXschemes for arbitrary convection-diffusion problems.In Chapter 5, the application of IMEX schemes to reaction-diffusion systems is con-sidered. A description of the two chemical reaction-diffusion problem is given and a linearanalysis of this problem is carried out to determine how it differs from the convection-diffusion problem. Based on this analysis and experimental evidence, efficient IMEXschemes are determined. This chapter concludes by describing how to choose IMEXschemes for reaction-diffusion problems.Conclusions are presented in Chapter 6. Suggestions as to the direction of futureresearch are also proposed.Chapter 2First and Second Order IMEX SchemesProperties of general linear multistep IMEX schemes for the system of ordinary differen-tial equationsit = f (u) vg(u)^ (2.3)are derived in the first four sections of this chapter. Subsequent sections determinestability properties of first and second order IMEX schemes for the one-dimensionallinear advection-diffusion equation. These results are compared to numerical experimentsof more complicated examples in Chapters 4 and 5.2.1 General Linear Multistep IMEX SchemesWe now derive s-step IMEX schemes for (2.3), s > 1. Letting k be the discretization stepsize and un denote the approximate solution at to = kn, these schemes may be written,11 8-1^s-i—kun-1-1 - Ea3 U^- E bif(un-i) v E cig(un-i)^(2.4).3=v^j=0^3=-1where c_ 1 0. For a smooth function u(t), expand (2.4) in a Taylor series about to = nAtto obtain the truncation error. This yieldss--1^s--1^kP-1^s-i1c [1 + Ea .ju(tn ) + [1 — Ejaja(tn ) +^+ 1^ [1 + E( —j) Pa3 1u (P) (t.)i=o j=i P.^j=is-1^s-i.^(I fdt^kp-1 S-1 dP-1^—E b3 f (u(tn )) + kE j b.; -`'2- It=tn^• (p 1).ID— i )P-11)3 dtP- 1 i t= t„—:7=0^i=1. J=18-1 s-i.^dg1—v E c 3 g(u(tn )) — kv[c-i — E c3]7,7 .t=t„j=-1 j=-18Chapter 2. First and Second Order IMEX Schemes^ 9^P-1^s-i^dP-ig^(p -1)!Ii[c_1 + E ci(—j)P-1]dtP-1 I,= tn + o(kP).^ (2.5)i=1.Applying Equation (2.3) to the truncation error (2.5), an order p scheme is obtainedprovideds_i1 + E ai^=^0i=os-i^s-i^s-i^1 -Ejai =^E b,^, E cij=i^i=o j=-11^s—i .2 s-1 s-12 + E j^—a 3. =^- E ibi , c_ i — E jca. ., 23=^j=1^ j=1(2.6)1 + S-1 iyai = S-1 ( jr1bi C-1^s1 E^P . j=i^j=i (p - 1)!^(p - 1)! jEr=i (p - 1)!Having determined constraints (2.6) for an order p, s-step IMEX scheme, we next examinehow to select p for a given s.2.2 Optimal Accuracy for IMEX SchemesWe begin by showing that the (2p+1) constraints of system (2.6) are linearly independentprovided p < s. Because linear independence for p = s implies linear independence forp < s, we need only consider the case p = s.The coefficient matrix for [a o , • • • as_1, bo, • • • bs-1, c-i, • • • cs_i]T representing the sys-tem (2.6) isAs^0. • • 00o (-1)' ...(1-s)'^-DsAs(2.7)0^-As^AsChapter 2. First and Second Order IMEX Schemes^ 10where A S and D s are s x s matrices,1^1^1^10^—1^—2^1 — s0^1^4^• • •^(1 — s) 20 (- 1) 8- 1 (-2)s- 1 • • • (1 — s)s- 112D s =The matrix A s is nonsingular. To see this, note that it is a Vandermonde matrix basedon the distinct values 0, —1, —2, ... , 1— s, and as such is nonsingular (see Golub and VanLoan[9]).Because A s is nonsingular, it is easy to see that columns 1, s+1,... , 2s, 28+2,^, 3s+1 of matrix (2.7) are linearly independent. Thus all (2p + 1) constraints are linearlyindependent and system (2.6) must admit a (3s — 2p) parameter family of solutions forp < s. Since column (23 + 1) of matrix (2.7) does not figure into the above proof weknow can have any nonzero value. This property verifies that the family of schemesis IMEX.We next verify that an s-step IMEX scheme cannot have accuracy greater than 0(k 3 ).Any 0(ks+r) scheme using s levels, r > 1, must satisfy the following constraints,s_iE bj = E ciJ-1E(—i)bj^+3=1^ j=1As =andChapter 2. First and Second Order IMEX Schemes^ 11(2.8)s1^ s-1E(—j)sbj =^E(—j)scij=1^ j=1Letting pi = cj bi the system (2.8) yields,s_ic_ 1 +^,aj = 0i=os--i+ E ,a3 ( —j) = 0j=1(2.9)s- 1c-i + = 0j=1Similar to the preceding proof the coefficient matrix for (2.9),1 1^1 11 0^— 1^• • — (s — 1)B1 0^(-1)s^• (1 — s)sis nonsingular and soC-1ItoB = 0its-1implies that c_ 1 = 0 and b.; = ci, j =1,...,s — 1. Since an IMEX scheme requires thatc_ 1 0, an s-step IMEX scheme is unable to achieve more than orders accuracy.Because it is desirable to achieve the greatest accuracy with the fewest steps, theremainder of this work only considers s-step, O(ks) schemes. As shown above, the familyof such schemes has exactly s free parameters.Chapter 2. First and Second Order IMEX Schemes^ 122.3 Method of AnalysisAll the IMEX methods (2.4) are consistent, provided the constraints (2.6) are satisfied.By the Lax-Richtmeyer Theorem (see Strikwerda[19]), these methods also are convergentwhen applied to a well-posed initial value problem, provided stability is ensured.General, sufficient stability criteria are difficult to specify because the ODE sys-tem (2.4) is potentially very large. To derive stability properties of IMEX methods,we consider the one-dimensional advection-diffusion equation, Ut = aUs vUsx where aand v are constants, v > 0, subject to periodic boundary conditions. Properties derivedfrom the analysis of this simple prototype problem should provide clues to the behaviourof discretizations of more complicated parabolic-hyperbolic systems.Using centred approximations, D 1 and D2, for the first and second derivatives, re-spectively, we obtain the corresponding semi-discrete equations,Ui = aD i Ui vD2 Ui, 1 < i < M.(Here, a uniform spatial grid with M points has been employed.) Applying a discreteFourier transform diagonalizes this system to= i/3 xi aixi, j = 1 . . . Mwhere cad and Oj are given in (2.13) below. For notational convenience, we write=^cex.^ (2.10)Applying the general multistep IMEX scheme (2.4) to (2.10) yields3-1^3-1^3-11^1xn+ 1 + — E aixn-3 = E ba ioxn-3 + E ciaxn-j (2.11)which is a linear difference equation with constant coefficients. (In comparing (2.11) to(2.4) note that v is buried in a.) The solutions of Equation (2.11) are of the formxn+1 =Pi Ti + P2'7T + • • • + PsT:Chapter 2. First and Second Order IMEX Schemes^ 13where rj is the jth root of the characteristic equation defined bys-i(I)(z) (1 — c_ i ak)zs +E( 3a — b i0k — ci ak)zs-j -1^(2.12)i=oand pi is constant for Tj simple and a polynomial in n otherwise. Clearly, stability holdsfor I Tj < 1, Ti simple and 1 731 < 1 where T3 is a multiple root.Throughout this thesis, we consider D 1 and D2 such that D i Ui = U'+'2handD 2 U, = tf,+1-2h2s+ut-1 This determines that2v^a(a.; , 03 ) = (— [cos(271-jh) — 11, —h sin(27rjh))^(2.13)which lie on the ellipse of Figure 2.1. Determination of time-stepping restrictions isaccomplished by finding the largest time step such that the ellipse lies in the absolutestability region of the IMEX scheme.Another important property can also be examined using these techniques. Fromthe exact solution of x = ii3x^ax it follows that lx(tn+i)I = eivIX(tn)1. Thus theChapter 2. First and Second Order IMEX Schemes^ 14amplification associated with the differential equation (2.10) over a time interval of lengthk is given by ea k . The corresponding contours in the a-/3 plane are plotted in Figure 2.2.From this figure, we note that, ideally, the roots of the characteristic polynomial for ourmethod should be small for a large and negative. This behavior corresponds to fastdecay of high frequency components when v is large. As discussed in Sections 1.2 and1.3, this behaviour often is particularly desirable when applying spectral collocation toan aliased computation or if a multigrid method is contemplated for the solution of theimplicit equations.2.4 Schur and von Neumann PolynomialsFrom the previous section we see that stability results for IMEX schemes are obtainedby examining the roots of the characteristic equation (2.12). To obtain results for the/3-axis for the higher order schemes considered in Chapter 3 we use the theory of Schurand von Neumann polynomials (see Miller[14]). The relevant definitions and theoremsChapter 2. First and Second Order IMEX Schemes^ 15from the discussion in Sections 4.2 and 4.3 of Strikwerda[19] are outlined below.Definition 2.4.1 The polynomial co is a Schur polynomial if all its roots, r j satisfyIT.ii < 1 .Definition 2.4.2 The polynomial^is a von Neumann polynomial if all its roots, rjsatisfy^1.Definition 2.4.3 The polynomial co is a simple von Neumann polynomial if co is a vonNeumann polynomial and its roots on the unit circle are simple roots.Definition 2.4.4 For any polynomial co , = Ri=o a izi the polynomial co* is defined byco* = >i_0 ds_iz i , where a denotes the complex conjugate of a.Definition 2.4.5 For any polynomial co o (z) we define recursively the polynomial (p j+i (z) =co (Op, (z)— co' (Op; (z)zTheorem 2.4.1 (pj is a simple von Neumann polynomial if and only if either Icoi(0)1 <1V;(0)1 and co3+1 is a simple von Neumann polynomial or co 3+1 is identically zero and co",is a Schur polynomial.Theorem 2.4.2 A necessary condition for the finite difference scheme to be stable is thatthe characteristic polynomial 0 defined in (2.12) is a simple von Neumann polynomial.2.5 First Order IMEX SchemesWe proceed by deriving one-step, first order IMEX schemes for Equation (2.3). Onerepresentation of this one-parameter family of schemes yields,un+1 un = k f (0) + vk[(1 7)g(un) ,yg(un+i )] —^ (2.14)Chapter 2. First and Second Order IMEX Schemes^ 16and we restrict 0 < -y < 1 to prevent large truncation error. As we have shown in Sec-tion 2.2, this one-parameter family must describe all first order, one-step IMEX schemes.Some of these schemes are familiar. For example, the choice 7 = 0 yields the ForwardEuler scheme,un+1 — un = k f (un) + I) kg(un).This scheme is fully explicit, and will not be considered further.Another possibility is to apply Backward Euler to g and Forward Euler to f. Thischoice (y = 1) yields,un+1 — un = k f (un) + v kg (un+ 1 ).^ (2.15)Backward Euler is the first order member of the class of backward differencing formulas(BDF's) (see Gear[8]). These implicit schemes are centred in time about time step (n+1).Implicit-explicit methods such as (2.15) which apply a BDF to g and which extrapolatef to time step (n + 1) will be referred to as semi-implicit BDF (SBDF) schemes.Having determined integration formulae, we next consider stability properties for firstorder IMEX schemes.2.6 Stability of First Order IMEX SchemesFrom Equation (2.14) the general first order IMEX scheme applied to x = (a + i /3)x isxn+ 1 = (a, 0)xnwhere (a, 0) = 1 + ka(1 — 7) + ik/31 — k-yaThe stability region is thus 1(a, 0) : 1(a, 0)1 5_ 1}.Figure 2.3 presents stability contours for first order SBDF, y = 1 Strong decay occursfor a large and negative. Furthermore, the scheme allows variable time-stepping and usesrelatively little storage. For first order schemes, the choice y = 1 is often preferred sinceChapter 2. First and Second Order IMEX Schemes^ 17smaller y values do not have as good decay properties and may become marginally stableor unstable for large v. See Figure 2.4 for the case 7 = 2 , which yields the second orderCrank-Nicolson scheme when f 0.All first order IMEX schemes have the shortcoming that they are unstable for v =0, a^0, sinceVO, = 11 -1- ik,31 = V1 + k 2 /3 2 > 1.Furthermore, at least a second order time-stepping scheme is often desirable since asecond order spatial discretization is used. We thus consider second order methods forthe remainder of this chapter.2.7 Second Order IMEX SchemesApproximating is = f(u) vg(u) to second order using IMEX schemes leaves two freeparameters. If we centre our schemes in time about time step (n -y) to second order,Chapter 2. First and Second Order IMEX Schemes^ 18we obtain the following family,_1{(,k), + _1 \ un+i - 27un + (,), _ 2_1 )un-i ] =^2 1^1(7 + 1) f (Un ) — -y f (2171-1 ) +^/1{(7 +^22— )g(un+1 ) + (1 - -y - c)g(un)-1- -c g(in-1)].^(2.16)Some of these methods are quite familiar. For example, selecting (y, c) = (li , 0) yields_I [un-1-1 - un] = _3 f(un ) _ _I f(un-i ) + _v [g(un+1) + g(un )] (2.17)k^2^2^2which is one of the more frequently used schemes in spectral methods applications (e.g.Canuto et al.[6] and Kim and Moin[13]). Because it applies Crank-Nicolson for theimplicit part and second order Adams-Bashforth for the explicit part, this scheme will bereferred to as CNAB (Crank-Nicolson, Adams-Bashforth). In the next section, we showthat the best asymptotic decay properties for y = 2occur when c = 8. This choice gives1_ run-Fi - unj = _3 f (un)^1 f(in-i) + v[ 9 g (un-Fi. ) + _3 g (un) + 1 g / un-].(^)]^(2.18)0^2^2 16^8^16Chapter 2. First and Second Order IMEX Schemes^ 19Because of the obvious similarity to CNAB, this scheme will be called modified CNAB.By setting (7, c) = (0, 1) in (2.16) we obtain another method which has been appliedto spectral applications (e.g. Brachet et al.[3]),2k [un+1 - unl = f (un) + 2 (un+1) + 9 (un-1 ) 1 •^(2.19)This scheme uses leap frog explicitly and something similar to Crank-Nicolson implicitly.For this reason, this method shall be referred to as CNLF (Crank-Nicolson, Leap Frog).Finally, setting (-y, c) = (1, 0) yields2k [Q un+ 1 - 4un + un -1 ] = 2f (un) - f (un -1 ) + vg(un+ 1 )^(2.20)which shall be referred to as SBDF since this scheme is centred about time step (n + 1).Other authors, such as Varah[22], call this scheme extrapolated Gear.Having determined integration formulae, we direct our attention to obtaining stabilityproperties for second order IMEX schemes.2.8 Choice of c for Second Order SchemesThe second degree characteristic polynomial resulting from (2.16) applied to x = (a+i )3)xis given by0(z) = fry + 2 _ ak(-y + -C )]z 21-[2-y -I- if3k(7 +1) + ak(1 - 7 - cAz + -y - -2 ii3k7 - ak-c-2(2.21)Because the parameter c in Equation (2.21) is always multiplied by a we choose c ac-cording to stability properties for lal >> 1/k. For this case, the roots of the characteristicequation (2.21) are given approximately by(7 + 2 )'T 2 + (1 - 7 - Or + 2 = 0Chapter 2. First and Second Order IMEX Schemes^ 20which gives usc —1 +(1— 7) 2 — 2cT1,22 7 CFor any (y, c), evaluatingDrys - max(Irii, ir21)^ (2.22)provides an estimate of the high frequency modal decay for large v. Minimization over cdetermines the method with the strongest asymptotic high frequency decay for a partic-ular 7.First consider c < (12)2 . Using elementary calculus, it is easy to verify thatV(1 — 7) 2 — 2cmin ^c< (1 —2-y)2^2y + candmin 1 7 + c —Ic< 0-2-,02 27 + coccur when c (1--y)2 Thus the minimum for /71„, must occur when c > (1--y) 2 -^2 • For2^•c > (1' )2 we have2 min 7311„c = minc> (1-02^ .(1-7)2 2 c' 21— y27+c^1+-y .(2.23)From which it is clear that the minimum is attained for c = (1 212 if 7 > 0 and c > 2 if7 = 0.In subsequent sections we use the following simple estimates of the high frequencydecay:1 — 7D7,0 =, 1 --y =Dry (1 —10221— y1+7'1— y+1 •Chapter 2. First and Second Order IMEX Schemes^ 21These results show that for c 0, we need -y^to ensure stability for all v. Furthermore,from Equation (2.23), Di„, = 0 <=> (7, c) = (1, 0) which implies that SBDF possesses thestrongest asymptotic decay of second order methods.2.9 Second Order Instability along the Nonzero ,3-AxisWe proceed by verifying stability at the origin, (a, ,3) = (0, 0). At the origin, the rootsof the characteristic equation (2.21) are 1 and 2.7-1 These are simple roots, whose2-y+1 4magnitude does not exceed 1 provided -y > 0. Thus by Theorem (2.4.2), (a„3)=(0,0) isstable for all second order IMEX schemes such that y > 0.Stability along the )3-axis is also a desirable characteristic. Schemes which exhibit thisbehaviour are better able to compute solutions for problems with very small v and forconvection problems with small artificial diffusion. For this reason, the remainder of thissection considers second order methods applied to the one-dimensional wave equation,= i/3x.Applying a second order IMEX scheme such that ey = 0 to the problem^i/3x yieldsthe leap frog scheme. Leap frog is known to be stable when applied to x = i/3x, providedk < —h7 as is outlined in many texts such as Strikwerda[19]. Thus the CNLF scheme, ina particular, is stable on the ,3-axis provided k < aFor other IMEX schemes, we know from Varah[22] that when a = 0 one of the roots,Ti, of Equation (2.21) satisfiesITi(#01 2 = 1 + (-y2 + )(0k) 4 + • • • > 1for /3k sufficiently small and -y > 0. Thus for all k, all second order schemes such thaty > 0 are unstable on the nonzero /3-axis.Chapter 2. First and Second Order IMEX Schemes^ 222.10 Second Order Stability ContoursFurther information can be obtained from the stability contours in the a - 0 plane. Theseplots are displayed in Figures 2.5 to 2.10. Figure 2.5 shows the contours for CNLF. Thismethod is stable for all v > 0, provided k < !. Such a time step restriction is undesirablesince it applies even for large v and small h. Furthermore, the decay of high frequencymodes can be weak, tending to 1 as a tends to —oo. Comparison of CNLF to othersecond order methods such that y = 0 suggests that CNLF produces the largest stabilityregion amoung such methods.The contours for CNAB are plotted in Figure 2.6. This method has a reasonabletime step restriction for larger v and small h. It is unstable near the /3-axis, however.It also suffers from poor decay of high frequency modes, since the decay tends to 1 as atends to —oo. Using modified CNAB, (7,c) = (,-1g.), the decay tends to 3 , a significantimprovement. See Figure 2.7 for these contours.Chapter 2. First and Second Order IMEX Schemes^ 23-1 0/k^ aFigure 2.6: Stability Contours for CNAB, (7, c) = (2, 0)- 1 0/k^ 0aFigure 2.7: Stability Contours for Modified CNAB, ('y, c) = (2, 8)Chapter 2. First and Second Order IMEX Schemes^ 24The contours for SBDF are displayed in Figure 2.8. This method has the mildesttime-stepping restriction when v is large and h is small. The decay of high frequencymodes is strong, tending to 0 as a tends to —oo. This method, however, has the strictesttime step limitation for small la l .For y = .25, stability contours for c = (17 )2 = .28125 and c = (1 —7) = .75 are plottedin Figures 2.9 and 2.10. Comparison of these contours indicates that c = (1 — y) = .75possesses a milder time step restriction for large v. Thus the method with the strongestasymptotic decay does not necessarily allow the largest stable time steps when v is large.Many other methods can be considered. The next section discusses which secondorder method allows the largest stable time steps for a particular problem.2.11 Choice of Second Order MethodWe now develop a quantitative method for describing time step restrictions for secondorder schemes. Such a method will help select which second order scheme to use for a.800n11111^it^II-10/k aFigure 2.10: Stability Contours for (-y, c) = (.25, .75)-11 1111111111111111111111111111111111^1111111111111111111111111^1111111111111111111111111111111• 0. 800.8001/kChapter 2. First and Second Order IMEX Schemes^ 25-10/k^ aFigure 2.9: Stability Contours for (-y,^= (.25, .28125)Chapter 2. First and Second Order IMEX Schemes^ 26particular problem.We begin by choosing a method (2.16) having fixed but arbitrary y and c. For theproblem is = aux + vuxx , a > 0, we apply this discretization scheme using spatial stepsize h and temporal step size k. Let the time step restriction on k be k. Similarly, forthe problem i) = v' v xx , a' > 0 define h', k' and k'. Suppose 4,,a = Jai. Then thescaled (2.13) gives(cVj' i@'.1) = a/a' /hh/ (h22v^a[cos(21-jh, ) —^T sin(27rj h'))which implies that the ellipse of Figure 2.1 is scaled by a factor 21,-g. Setting k' =^kcauses a similar scaling in the stability contours. Thus /22' — ci and so it is consistenthi^hto plot ha ^a function ofThe quantity hrepresents the ratio of viscous effects to convective effects betweenmesh points. The inverse quantityahR = (2.24)has been referred to as the mesh Reynolds number (see Peyret and Taylor[17]).As can be seen from Figures 2.11 to 2.13, increasing y allows larger stable time stepswhen h> 2. The case c = 1 — y also has the property that decreasing y allows largertime steps for 1-cfr, < 2. Comparison of Figures 2.11 to 2.13 indicates that the largesttime step can be applied using SBDF for h> aand CNLF for h< 2. This resultphysically corresponds to selecting SBDF when discrete diffusion dominates, and CNLFwhen discrete convection dominates. From this perspective, the popular CNAB is onlycompetitive when vFrequently, an important consideration when choosing a second order scheme is whatthe constant of the truncation error is. For example, Crank-Nicolson is known to have amuch smaller truncation error than second order BDF (see Hairer et al.[11]), so we expectCNAB to have a smaller truncation error than SBDF. Modified CNAB is expected to10 -211.1111^ I^I^.11111 10-1^100^ 101Vracosity/ahFigure 2.11: Time Step Restriction For Various -y where c = 1 — -y10 11` 0a10 -1Chapter 2. First and Second Order IMEX Schemes^ 27have a truncation error similar to CNAB, however. (Numerical experiments in Chapter 4support this claim.) Because of its small truncation error and because it produces strongerhigh frequency spatial decay than CNAB, modified CNAB may be preferred in certainproblems over CNAB or SBDF when h> 2. To obtain a better understanding of therelative efficiencies of second order IMEX schemes, a study of the truncation errors wouldbe useful.Although the results of this section can be applied to a wide variety of problems, westill seek a method which is stable for all v > 0 and has strong decay for Ic >> 1/k. Toachieve these objectives the next chapter examines higher order IMEX schemes.Chapter 2. First and Second Order IMEX Schemes^ 281 0210-2^ 10-1^ 100^10 1Viscosity / ahFigure 2.13: Time Step Restriction For Various 7 where c = 010 1Ls 1 0 01 0 -110-2^10-1^ 100Viscosit y /ahFigure 2.12: Time Step Restriction For Various -y where c = (1 7 )21010 11 0 110 -1..N i o°4JeChapter 3Higher Order IMEX SchemesNone of the second order IMEX schemes have all the stability characteristics we desire.CNLF is stable for all v > 0 in the linear case, but has poor decay for large lal. SBDFhas good stability characteristics for a << -1/k but has a very strict time step restrictionfor small lal. Furthermore, no two-step, second order method contains the 13-axis withinits stability region.To find a scheme satisfying all the desired stability properties we now consider higherorder IMEX methods.3.1 Third Order IMEX SchemesWe begin by deriving third order, 3-step IMEX schemes for it = f (u) + vg (u). FromSection 2.2, these schemes form a three parameter family of methods. One possibleparameterization yields,1 ^ 3[( _ 72^+ 7 + 1_ +09)un+i + (__72 - 27 + 1- - 19)un +k 2 3^2^23 1^1( _272 + 7 _ nun-1 + (__272 + _d un-2] =( 72 + 31' + 1 + 23 0)./(0) (72 + 27 + 4 6)f(un -1 ) + (72 + 7 + fle)f(un- 2 ) +2^12^ 3^2^1214( 7 2 + 7 + og(un+i ) + (1 112 3c + 2320)g ( un) +2 1^'^,v 22^,,,^4-3( 7^ ' + 3c - 61)g(u'^12l) + (-5 61 - c)g(un- 2 )]. (3.25)29Chapter 3. Higher Order IMEX Schemes^ 30These schemes are centred about time step (n +-y) to third order, provided 0 = 0. As forlower order schemes, the value of -y should be between 0 and 1 to avoid large truncationerror. Also, the parameter c is multiplied by v, so this parameter should be adjusted tomodify large v properties of a scheme.Some of these methods are familiar. Letting 0 -p +oo yields23 ^n-2)](un+1 un) .= {f (un) g(un )] _3 [f(un-i ) g(un-i )1 +12 [f(un- 2) g(uwhich is the third order Adams-Bashforth method.Setting (-y, 0, c) = (1, 0, 0) yieldsk 6^2^3^) (3.26)which is the third order BDF for the implicit part, and similarly to Section 2.5 is calledthird order SBDF.Having determined integration formulae, we direct our attention to obtaining stabilityproperties for third order IMEX schemes.3.2 Choice of parametersThe third degree characteristic polynomial resulting from Equation (3.25) applied to= (a + it3)x is given by24)(z) ..=[272 + 7 + 31 + _0 ( 7 2+ 7 + c)alciz32- [1y 2 + 27 - 1 + 0 + ( 7 + 37 + 1 + 23 0)i/3k (1 7 2^230 k] 22^2^2^U 1 + - 7 - 3c + —12 )ct z_ 72^4+ [-37 2 + 7 - 1 + (-y 2 + 27 + -40)i i3k + (7 2^3c + -3 0)aldz2 3 2- [12 - 6-1 + (7 + 7 + 2^12—5 0)iOk + (-5 0 - Oak]. (3.27)2^_1 ( 11 un-Fi 3un 3 un-1 _ 1un-2) = 3f(un) - 3f(un-1 ) Aun-2) vg(un+iWe now determine which methods produce the strongest asymptotic decay as a -oo.For this case, the roots of the characteristic polynomial (3.27) are given approximatelyChapter 3. Higher Order IMEX Schemes^ 31by2 + 7( 7 " + C) Z3 + (1 72 3c+ 23 0)z 22^ 127 y21^3c + -4O)z + 50 - c = 0. (3.28)2 3^12By determining the solutions, ri , 72 and 'r3 , of Equation (3.28) we may evaluateRy , e ,,^maX ( 1 71 I , 1 72 1 1 1 731)to obtain an estimate of the high frequency model decay for large v. Minimizationover (7, 0, c) determines the method with the strongest asymptotic high frequency decay.Certainly if1)10,00, C° = 0^(3.29)then (70 , 00 , co ) minimizes the amplification as a --4 -oo. From Equation (3.28) we satisfy(3.29) ifEach term is divided by ( 1 ---12-1---M2by letting ( -Y2-Q 1—' + co ) -* ±00.1 - I'd - 3c0 + E 00 = 0-ro + A)- + CO2^2C)- + 3c0 — 1°0 = 0^ (3.30)-Y1, -Flo 2 + COMO — CO2^ = 0./0-10 2 + CO+ co ) to allow the possibility of satisfying Equation (3.29)Simplification of the system (3.30) yields(3 -Yo - 1)(yo - 1)=02411 d- co00 — 6(7(1 -- 70) -Y02 +-V° = 02 + COn172- v CI — c0 = O.„2_Lr 0 -CY 0 + CO(3.31 )(3.32)(3.33)Chapter 3. Higher Order IMEX Schemes^ 32For C2--Q4L° -I- co ) finite there are two possibilities, (yo , 00 , co) = (1, 0, 0) and (yo , 0o , co) =( i , _1 , __:), both of which specify third order SBDF. For ( 2L-42 —' + co ) infinite, Equa-tion (3.31) implies 'To << 'col. Using this in Equation (3.32) implies 100 1 << 'co l. Applyingthese results to Equation (3.33) results in a contradiction, so C 2-±c' + co ) must be finite.Because third order SBDF has the strongest asymptotic decay of third order IMEXschemes, special attention is given to its properties throughout the remainder of thischapter.3.3 Third Order Stability along the /3-AxisWe begin by verifying stability at the origin (a, )3) = (0, 0). From Equation (3.27) theroots of the characteristic equation at the origin are given by1^1^5^ 1( T - 1 )[(7 2 + 7 + i + 61 ^+ (7 2 + 7 — —6 )T + .1 -Y2 — —6 ] °To be stable all roots must have magnitude less than or equal to 1, and those havingabsolute value 1 must be simple. Third order SBDF, as well as the methods of Section 3.4satisfy this criteria.Stability along the /3-axis is verified by considering the application of third order meth-ods to the one-dimensional wave equation. We demonstrate this technique for third orderSBDF. Applying Equation (3.26) to x = i/3x we obtain the characteristic polynomial3^1(p0 (z) = —16 z3 — 3(1 + i/3k)z 2 + (-2 + 3i/3k)z — —3 — i/3k^(3.34)Using the definitions of Section 2.4 we determine thatVo'(0)1 2 = 3.3611 > .1111 + 1.0000/3 2 = koo (0)1 2for 1/31 sufficiently small. Using the symbolic computation package Maple[7] we find thatSo l (z) 0 0 and that144(0)12 = 10.5625 — 6.5000/3 2 + 1.000004 > 3.0625 + 1.7450/3 2 + 9.00000 4 = 1',91(0)12Chapter 3. Higher Order IMEX Schemes^ 33for 131 sufficiently small. Reapplying Definition 2.4.4 we determine that the square of theabsolute value of the root of cp 2 (z) is given by56.25 — 123.75/32 — 51.9434 + 132.00/3 6 + 64.00 /3856.25 — 123.75 /3 2 — 88.5034 + 234.25/3 6 + 36.0038which is less than 1 on the nonzero /3-axis near the origin. Thus yoo (z) is a von Neumannpolynomial by repeated application of Theorem 2.4.1, and so Theorem 2.4.2 assures usthat third order SBDF is stable on the /3-axis for sufficiently small.Similar studies for the methods described in the next section indicate that they alsoinclude the /3-axis in the absolute stability region.3.4 Third Order Stability ContoursFurther information about stability in the a-3 plane can be obtained by plotting max{ Izi :0(z) = 0}, where 0(z) is defined in Equation (3.27). These stability contours are dis-played in Figures 3.14 to 3.19.We begin by examining if it is possible to arrive at a stable scheme for any fixed -y.For a fixed, but arbitrary y and for 1 0 1, ci oo we obtain an approximate local minimumof D.1„ 0,, if e = 0.4661. Using these parameters, the scheme simplifies to_1 (un+i — un) 23^4^512 f (un) — f (un —1 ) + 2- f (u') + .4661vg(un+ 4 )+.5184vg(un) .0650vg(un -1 ) — .0494vg(un -2 ) (3.35)which applies third order Adams-Bashforth to the explicit term. The stability contoursof Figure 3.14 suggest that this method is stable for all v > 0 provided k < 0.62 Thisrestriction is more severe than that for third order Adams-Bashforth applied to x = i,3xbecause of the dip in the stability contours when a < 0. Careful analysis similar to thatof the previous section verifies that the 3-axis is included in the absolute stability region.Chapter 3. Higher Order IMEX Schemes^ 34This result demonstrates that for third order methods, it is possible to find methods forany -y which are stable for all v > 0 by varying 0 and c.In Chapter 2, the most interesting second order schemes were produced by selecting7 equal to 0, 2 or 1. We consider schemes for these values of -y below. The parameters 0and c are chosen to produce schemes which allow large stable time steps as I) —) oo.For example, the method (7, 0, c) (0, —2.036, —.876) of Figure 3.15 is stable for allv > 0 provided k < 0.67! . Similarly (7, 0, c) = (.5, —1.21, —.5) of Figure 3.17 is stablefor all v > 0 if k < 0.65!. In both these cases substantially larger time steps can betaken for large or moderate jal than for the method (3.35). Furthermore, stronger highfrequency decay occurs for these methods than for method (3.35).From Section 3.2 we know that the strongest decay as a —> —oo occurs for thirdorder SBDF. The stability contours for this method are shown in Figure 3.18. This plottogether with the zoom-in of Figure 3.19 suggest that third order SBDF is stable for alla < 0 provided k < 0.62L'a . The plot of Figure 3.18 also indicates that very large timesteps can be taken for large or moderate Although applying -y = 1 and nonzero 0and c may allow somewhat larger stable time steps we focus on 0 c 0 since otherchoices degrade the strong asymptotic decay.For -y 1, the choice 0 = c 0 is not recommended. As seen in Figure 3.16,(y, 0, c) 0, 0) results in a small stability region. This plot indicates that very smalltime steps are needed for moderate or large since the ellipse of Figure 2.1 must lie inthe stable region.'kChapter 3. Higher Order IMEX Schemes^ 35Chapter 3. Higher Order IMEX Schemes^ 36a41 1111111 11111111 11 111111 111111111111111111111111111111111111I ^11111111-9/k=- 1 .00 1. 00900^9,0G1 a41111^111111111111111111^111111111111N111111111111111111111111111111111111111111111111111111.1-ea800,-1/k• 600. 600-= . 400I I 1 1 111ii 0Figure 3.18: Stability Contours for Third Order SBDF, (7,0, c) = (1, 0,0)4-1.11111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111-==1/k-.2/k^4111111111111111111111111111111111111111111111,1111111111111111111 ^111111111111111111111111111111117^0^0 .1/kaFigure 3.19: Zoom-in Around a = 0 for Third Order SBDFChapter 3. Higher Order IMEX Schemes^ 37Chapter 3. Higher Order IMEX Schemes^ 383.5 Fourth Order IMEX SchemesFrom Section 2.2, the general fourth order, four-step IMEX scheme is a four parameterfamily of methods. From the previous section we anticipate that the fourth order semi-implicit backwards differencing formula (fourth order SBDF) may have good stabilityproperties. For it = f (u) vg(u), this scheme is given by1 25(un+ 1 - 4un + 30 -1 - 4—0-2 1+ — Un-3 ) =k 12^3^44f (un) - 6f (un -1 ) + 4f (0 -2 ) - f (u"-3 ) vg(u"+1 ).^(3.36)The characteristic polynomial obtained from applying Equation (3.36) to x = (a + i/3)xis1^(z) = (-2512^ 3- ak)z4 - (4 + 40k)z3 + (3 + 6if3k)z2 - (-4 + 401v)z + -4 + i# k.Clearly, all of the characteristic roots tend to 0 as a -3 -oo so the amplification musttend to 0 as a -> -oo. Similar to the proof of Section 3.3, repeated application ofTheorem 2.41 using Maple[7] shows that the /3-axis is included in the stability region.Furthermore, the characteristic stability contours plotted in Figure 3.20 indicate thatfourth order SBDF is stable for k < 0.52-h= and that larger time steps are permitted as vaincreases.Comparison of Figures 3.18 and 3.20, however, indicates that the stability region issmaller than for the third order case, so smaller time steps may be required. Furthermore,comparison of Figures 3.19 and 3.21 shows that the /3-axis is closer to the boundary of thestability region for fourth order SBDF, suggesting that third order SBDF may dissipateerror better for v ti 0.Chapter 3. Higher Order IMEX Schemes^ 39Chapter 3. Higher Order IMEX Schemes^ 403.6 Comparison to Lower Order MethodsThird and fourth order SBDF methods are good choices for IMEX schemes for someproblems. For all ii > 0, these methods are stable for a time step restriction similarto the CFL condition. Greater accuracy and strong high frequency decay also makethese methods very attractive. Nonetheless, for many problems second order methodsare preferred. Higher order methods require more storage, and more work per time step.This extra expense could be justified if greater accuracy permitted larger time steps.Third and fourth order schemes, however, have more severe time step restrictions thansecond order schemes. In fact, Figure 3.22, shows that larger stable time steps can betaken with second order SBDF when -fr, > .11 for the linear advection-diffusion problemusing second order spatial discretizations. CNLF allows larger stable time steps thaneither third or fourth order SBDF forah < 1.Third order SBDF should be efficient for problems which require strong decay forlal >> 1/k and a moderate time step restriction for c''7, < 0.1. It should also be effectivefor problems where -fr << 1, since the /3-axis is within the stability region.Fourth order SBDF has a particularly severe time step restriction for the advection-diffusion problem when -,fri is moderate or large. For example, when ;Li = 8, fourth orderSBDF only can apply one tenth the time step of second order SBDF, as can be seen fromFigure 3.22. This restriction on the step size would appear to limit fourth order SBDF toproblems where accuracy, and not stability is the reason for limiting the time step size.Chapter 3. Higher Order IMEX Schemes^ 41Chapter 4Convection-Diffusion ExperimentsThe previous chapters have dealt with stability properties of IMEX schemes for theone-dimensional linear constant coefficient advection-diffusion equation. These resultsprovide necessary, but not sufficient conditions for stability for variable coefficient andnonlinear convection-diffusion problems.The next sections summarize numerical experiments which verify our analysis forthe simple, linear problem can be useful for determining which IMEX scheme to applyto more complicated problems. In particular, strongly damping schemes are shown tobe more efficient in certain spectral collocation and multigrid applications. In order tocalculate starting values for multistep IMEX schemes, we use one-step (low order) IMEXschemes with a very small time step, unless otherwise noted.4.1 Finite Difference Approximations4.1.1 Variable Coefficient ProblemTo examine nonzero viscosity behaviour, we consider the one-dimensional variable coef-ficient problemut sin(271- x)ux = vuss^(4.37)subject to periodic boundary conditions on the interval [0,1] and initial conditionu(x, 0) = sin(2irx).42Chapter 4. Convection-Diffusion Experiments^ 43Because shocks or other discontinuities do not form for any v > 0 for this problem, weexpect that centred differences can be used to compute accurate approximate solutionsfor this example provided a sufficiently fine spatial mesh is applied. For this reason,the next two sections use second order centred differences D 1 and D2, as defined inSection 2.3, to approximate us and uss .4.1.2 Small Mesh Reynolds Number ComputationsTo test the theory's predictions for small mesh Reynolds numbers (2.24), the modelproblem (4.37) was approximately solved using discretization step sizes h = a- in spaceand k = 1.8h in time. Use of such step sizes is appropriate only for strongly dampedflow. Utilizing several IMEX schemes, computations to time t = 2 are performed forviscosities, v, in the range .01 < v < .1. These values correspond to mesh Reynoldsnumbers, R, in the range 1.59 > R > 0.159.From Figures 2.11, 2.12, 2.13 and 3.22 the theory predicts that for these step sizesSBDF is stable for a larger viscosity interval than any other scheme. As v decreases, thirdorder SBDF followed by CNAB should become unstable. Stability for modified CNABshould be similar to, and somewhat better than, CNAB. Furthermore, fourth order SBDFand CNLF should be unstable over the entire interval because the theoretical stabilityrestriction is violated.By comparing results to those for h = *1 and k = .225h using SBDF, the max normrelative errors for second order schemes are evaluated. Third and fourth order schemesare compared to computations with these same step sizes but using third and fourthorder SBDF. The resultant errors are plotted against v in Figure 4.23. Fourth orderSBDF and CNLF are not included because they are unstable over the entire interval.The plots of the figure clearly coincide with the results of the theory.Figure 4.23 also indicates that SBDF is the only stable method for the above choiceChapter 4. Convection-Diffusion Experiments^ 44of discretization parameters when 0.0015 < v < 0.0025. This agrees with the predictionthat SBDF allows the largest stable time steps for small mesh Reynolds numbers (seeSection 2.11). Although third order SBDF has a smaller stability interval, it may beuseful in problems where high accuracy is needed since it produces a smaller error thansecond order methods when stable.For the second order family of schemes (2.16) such that ey = 4, plots for several c aremade. These are displayed in Figure 4.24. Although the best asymptotic decay occursfor c = .28125 (see Section 2.7), a wider variety of viscosities can be computed with otherchoices of c. An interesting and open problem would be to determine which second orderscheme for a given y permits the largest stable time steps when v >> 1/k. The choicec = 1 — -y appears to allow particularly large time steps for v 1/k.To conclude, Example (4.37) supports several of the theoretical results for small meshReynolds numbers. In particular, SBDF seems to be a good choice when the largeststable steps are desired and third order SBDF appears appropriate when higher accuracyis needed. Results for the large mesh Reynolds number case are considered next.4.1.3 Large Mesh Reynolds Number ComputationsTo test the theory's predictions for large mesh Reynolds numbers, example (4.37) wassolved using discretization step sizes h = 81 in space and k = .9h in time. Using severalIMEX schemes, computations to time t = 2 are performed for viscosities, v, in the range.001 < v < .01. These values correspond to mesh Reynolds numbers, R, in the range12.3 > R > 1.23.From Figures 2.11, 2.12, 2.13 and 3.22 the theory predicts that for these step sizesonly CNLF is stable over the entire viscosity interval. As v decreases, third order SBDFfollowed by SBDF and finally CNAB should become unstable. Stability for modifiedCNAB should be similar to CNAB. Furthermore, fourth order SBDF should be unstableChapter 4. Convection-Diffusion Experiments^ 45Chapter 4. Convection-Diffusion Experiments^ 46over the entire interval because the theoretical stability restriction is violated when v <.05 for these step sizes.By comparing results to those for h = -641 8 and k = .225h using SBDF, the max normrelative errors for second order schemes are evaluated for these computations. Third andfourth order schemes are compared to computations with these same step sizes but usingthird and fourth order SBDF. The resultant errors are plotted against v in Figure 4.25.Fourth order SBDF is not plotted because it is unstable over the entire interval. Theplots of the figure support the results of the theory.Figure 4.25 also indicates that CNLF is the only stable method for v < 0.002. Thisagrees with the prediction that CNLF allows the largest stable time steps for large meshReynolds numbers (see Section 2.11). CNLF is a particularly attractive choice becauseit has the smallest error of second order methods. Use of SBDF is not recommendedbecause it has the smallest stable interval and the largest error of any second orderscheme.For the same problem using k = .5h, dramatically different results are predicted be-cause all the methods satisfy their theoretical stability restrictions. Indeed, computationsfor CNLF, third and fourth order SBDF all produce errors which nearly coincide, sincespatial error dominates the solutions. Other second order methods appear stable andproduce only slightly less accurate results.To obtain more information about the behaviour of the solution for v near 0, we nextexamine a problem for which v = 0.4.1.4 Zero Viscosity ComputationsTo examine the limit v 0, we consider the one-dimensional nonlinear problem1u t + —2 cos(27rt)(1 u)u, = 0 (4.38)Chapter 4. Convection-Diffusion Experiments^ 47having periodic boundary conditions on the interval [0,1] and initial conditionsu(x, 0) = sin(27rx).As in the previous two sections, we use second order centred differences to approximateus and uss . For h = so and k = .5h, computations are performed to time t = 100. Thetime step value k = .5h was used to ensure that fourth order SBDF satisfied the stabilityrestriction k < .52h. Using the fact that the exact solution to this problem at integer tequals the initial data, i.e.u(x, n) = sin(27rx), n = 0, 1, 2, . . .we compute the error in the solution at time t = n, n = 1, 2, ... , 100. The relative maxnorm errors for several schemes are plotted in Figures 4.26 and 4.27.All second order schemes tested, such that y > 0, were unstable. This result supportsthe conclusion of Section 2.8. CNLF, third order SBDF and fourth order SBDF were allstable.Chapter 4. Convection-Diffusion Experiments^ 48The error at t = 100 for CNLF was less than 10'. This method performs extremelywell for this example since it is a time reversible method. It is known that leap frog isonly stable for v = 0, however. Thus CNLF may be more sensitive to perturbations ofthe eigenvalues away from the /3-axis than third or fourth order SBDF.As can be seen in Figure 4.27, the error for third order SBDF is concave down as tincreases and has a final value of 4% at t = 100. The fact that the error is concave downsupports the conclusion that the /3-axis is included in the stable region, since error isbeing dissipated.For fourth order SBDF, the theory predicts that the /3-axis is stable, but that itis much closer to the boundary of the stability region. Examination of the error tofour significant digits verifies that dissipation occurs for fourth order SBDF. Becausethe dissipation is very weak, the growth of the error appears linear as can be seen inFigure 4.27.In summary, CNLF, third and fourth order SBDF perform well for the zero viscositycase. Of these schemes, CNLF is non-dissipative while third order SBDF is the moststrongly dissipative. To obtain more information about the zero viscosity case, a linearanalysis of the phase lag error would be a useful undertaking.The next few sections consider problems with a stronger physical motivation.4.2 Spectral Collocation4.2.1 The Burgers EquationFor our next application we consider the Burgers equation,u t uus = vuss^ (4.39).0012 L,03O_c.0008U-L.0006L a0.0004cc.000209▪ .030Oo• . 025_cI— .0200Lu .0150•_•••.• 01 0.0050.50.45.40.35(ci .30LLd?. .25ct, .20.15. 1 0.05 5^10^15^20^25^30^35^40Figure 4.26: Zero Viscosity Example For Various Second Order Methods, y.040.0014.035Chapter 4. Convection-Diffusion Experiments^ 490^10^20^30^40^50^60^70^80^90^100Figure 4.27: Zero Viscosity Example For Higher Order MethodsChapter 4. Convection-Diffusion Experiments^ 501.0.8.6.4.27^0-.2-.8- 1 .0-1.0^-.8^-.6^ -.2^0^.2^.4^.6^.8^1 0Figure 4.28: Burgers Equation for Various t for v = 0.01subject to periodic boundary conditions on the interval [-1,1] and initial conditionsu(x, 0)^sin(irx).We expect that spectral solutions for this equation provide information about numericalsolutions to more complicated problems such as the Navier-Stokes equations.A plot of the solution of the Burgers equation for v = .01 at several different timesis provided in Figure 4.28. This computation was produced using Chebyshev collocationwith 40 basis functions and k = 1/160 using SBDF.The next few sections discuss Fourier and Chebyshev collocation implementations forthe above model problem. See Canuto et al.[6] or Boyd[2] for details on these methods.Chapter 4. Convection-Diffusion Experiments^ 514.2.2 Fourier Spectral CollocationSince the problem of this section is periodic, we expect that Fourier series makes a goodbasis of trial functions for this problem. Indeed, since1 au2^a2 UUt =V2 ax^ax 2we see that u t is antisymmetric for u antisymmetric. By selecting an initial conditionwhich is antisymmetric we guarantee that u remains antisymmetric for all t. Since onlythese components of the series contribute to the solution we use a Fourier sine series. Wethus approximate u by the seriesNUN(X,t) = E ai(t)sin(jx).j=1To determine ai (t), we enforce the differential equation at the collocation points,[auN^auN^a2uNiat +UN^ = V^ax^aX2 X= Xiwhere cr i (0) = 1 and ai (0) = 0, j^1. This scheme is also called a pseudospectralmethod since the nonlinear convection term is evaluated in physical space.By applying Fourier sine collocation with 40 basis functions and k = 1/40, we solvethe model problem at each time step to time t = 2. Because the system is small, theimplicit equations are solved in physical space using LU decomposition. In larger systemswhere greater efficiency is needed these would be solved in Fourier space using transformmethods (see Canuto et al.[6]). For CNAB, modified CNAB, SBDF and third orderSBDF the max norm relative error is plotted against viscosity (see Figure 4.29). Forsecond order schemes, the "exact solution" is based on a computation using N = 80modes and k = 3200 3200using SBDF. A computation using N = 80, k = and third orderSBDF provides the "exact solution" in the third order case.Chapter 4. Convection-Diffusion Experiments^ 521 0010- 1le-2O-I 10-310 -410 -5 10 -3^10 -1vi..0.ityFigure 4.29: Fourier Spectral Collocation for Burgers EquationCNLF was not included because the theoretical stability restriction is violated. Thiscan be easily seen because the linear advection-diffusion equation has eigenvalues (ian7r —vn2 7 2 ) for eigenfunctions e''. From these eigenvalues we know that the stability re-striction is k < —1 which is violated initially because lu(x, 0)1 00 = 1.As expected, third order SBDF has the smallest error of any of the methods whenstable. The stable region, however, is smaller than that for CNAB or SBDF. CNAB andmodified CNAB are once again very similar in behaviour with the modified version beingmarginally better. SBDF appears to allow the largest stable time step when v ti 0.01.Further refinement of the time step to k = —161 0 leaves third order SBDF as the methodof choice over the entire interval. Such a refinement seems unnecessary in this examplebecause the error is less than 1% for a step size which is 4 times larger.Chapter 4. Convection-Diffusion Experiments^ 534.2.3 Chebyshev Spectral CollocationBecause 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 (4.39) subject to the homogeneous Dirichlet bound-ary conditions u(+1, t) = 0, using a pseudospectral Chebyshev collocation scheme. TheGauss-Chebyshev grid,x i^cos[(2i — 1)r/2N], i = 1,^, Nis used for collocation points.Similar to the Fourier case, the max norm relative errors are evaluated for severalIMEX schemes using k = 1/40 and N = 40. As can be seen from Figure 4.30, the resultsare qualitatively similar to those of the Fourier case. SBDF performs particularly wellfor smaller viscosities. From the theory, it has the widest stability region for largeand thus is best able to accommodate the rapidly growing eigenvalues associated withChebyshev collocation.Both Chebyshev and Fourier collocation can be affected by aliasing (see Canuto etal.[6]). We consider aliasing for the Fourier case next.4.2.4 Aliasing in Pseudospectral MethodsAliasing occurs when nonlinear terms produce frequencies that cannot be represented inthe basis, and thus contribute erroneously to lower frequencies.For example, in the Burgers equation, the Fourier sine mode sin(mx), when explicitlyevaluated in the convection term produces the contribution sin(mx) 8sinatnx) = 2 sin(2mx).If 2m is greater than the number of Fourier sine basis functions, N, this frequencycannot be represented correctly and aliasing occurs. We now proceed to demonstratethat this behaviour can plague poorly spatially resolved computations when applyingweakly damping IMEX schemes.la lLO10SBDFa10 010 -1Ole -210 -40.SBDFI 10 -2Viscosit yur510 -3 10-1I^I^IChapter 4. Convection-Diffusion Experiments^ 54Figure 4.30: Chebyshev Spectral Collocation for Burgers EquationWe compute solutions for the model Burgers equation (4.39) subject to periodicboundary conditions and the initial conditionsu(x, 0) = 0.98 sin(27x) HF(x)whereHF(x)^0.01 sin(617rx) + 0.01 sin(627rx).We use Fourier sine collocation as described in Section 4.2.1 with N = 64 basis functions,and integrate to time t = 2.The value of the approximate solution at time t = k is obtained using first orderSBDF with the same step size. To represent the type of high frequency informationthat 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 firstorder SBDF step. Subsequent steps are then computed using the relevant second orderModified CNAB 3rd order SBDFk CNLF CNAB SBDFChapter 4. Convection-Diffusion Experiments^ 55scheme. For third order SBDF, the value at t = 2k is also needed. For the purposeof demonstrating aliasing effects, this value is computed using CNAB since we wish toretain most of the added high frequencies.The max norm relative errors for several IMEX schemes are computed by comparingresults to those for SBDF using N = 128 and k = —241N• In the third order case, resultsare compared to those for third order SBDF using N = 128 and k = —241N• These errorsare summarized below.1-64 •92 .575 .0056 .0072 .0045192 •060 .022 .0013 .00079 .0010For the case k = 4 we note that the error for CNAB is far greater than for modifiedCNAB, SBDF or third order SBDF. Using CNAB, a non-aliased computation usingthe 2/3's rule as described in Section 1.3 and k = s4 in a relative error of lessthat 10'. This non-aliased result, along with its aliased counterpart, are plotted inFigure 4.31. The "error" curve in this figure is that of the aliased CNAB. From thefigure it is clear that the main component of the error is proportional to sin(7rx). Thislow frequency mode is not part of the exact solution, and must be an aliasing artifactfrom high frequency components. It is interesting to note that even after 128 time stepsthe numerical solution is still plagued by high frequency components which have not yetdecayed. A similar study of CNLF reveals that it suffers from aliasing as well.Resorting to a smaller time step, k = 192, makes a substantial improvement in thesolution for CNAB and CNLF. Even so, the aliasing error for CNLF is sufficiently largethat further refinement is likely required.We conclude that use of a strongly damping scheme such as SBDF, modified CNABChapter 4. Convection-Diffusion Experiments^ 56.00050.00040.00030.00020.000100-.00010-.00020-.00030-.00040-.00050-1.0^-.e^ 0^.2^.4^.6^.8^1 0Figure 4.31: Aliased and Non-Aliased Computations for CNABor third order SBDF gives an inexpensive method to reduce aliasing in poorly resolvedcomputations. Application of weakly damping schemes like CNAB and especially CNLFmay necessitate undesirably small time steps to produce the high frequency decay neededto prevent aliasing in an aliased computation. Alternatively, weakly damping schemesmay be used in conjunction with an anti-aliasing technique such as the 3/2's rule orthe 2/3's rule. However, as described in Section 1.3, these anti-aliasing techniques ei-ther increase the expense of the computation or produce a severe loss of high frequencyinformation at each time step.As we shall demonstrate in subsequent sections, strongly damping schemes can be de-sirable for solving the implicit equations in finite difference equations when using multi-grid.Chapter 4. Convection-Diffusion Experiments^ 574.3 Time-Dependent Multigrid4.3.1 The Convection-Diffusion ProblemThe next model that we consider is the 2D convection-diffusion problem,ut (u • V)u = vAu^ (4.40)where u^(u, v). This model incorporates some of the^major ingredients of the 2DNavier-Stokes equations. Indeed, it is often being solved as part of projection schemesfor incompressible flows (see Hirsch[12]). We carry out our computations on the squaref2_-[0,1]x [0,1] and consider periodic boundary conditions and initial conditions,u(x,y,0)^sin[2r(x y)]v(x, y, 0) = sin[27r(x^y)].For spatial discretization we use standard centred differences. Thus we use A i„ thestandard 5-point Laplacian, to approximate A andwhereVh ==Dy^=DD:— ui_ i2h—2hto approximate V. For IMEX schemes the convection term, (u • V)u, is handled explicitlyand the diffusion term, vAu, is handled implicitly.This treatment yields a positive definite, symmetric, sparse linear system to solveat each time step. Such systems are solved efficiently using a multigrid algorithm, thecomponents of which are outlined in the next section.Chapter 4. Convection-Diffusion Experiments^ 584.3.2 Components of the Multigrid AlgorithmTo solve the implicit equations we apply a Full Multigrid (FMG), Full ApproximationScheme (FAS) algorithm (see Brandt[4]). Rather than applying the standard algorithmat each time step, we use the recent ideas of Brandt and Greenwald[5]. The main com-ponents of the resulting algorithm are outlined below.Smoothing is accomplished using Red-Black Gauss-Seidel. This relaxation techniqueis chosen because it has a very good smoothing rate. Prolongation is accomplished usingbilinear interpolation, and restriction by full weighting.The standard FMG cycle is modified so that the first coarse grid correction is per-formed before any fine grid relaxation. Based on the assumption that the incrementbetween time steps is smooth, Brandt and Greenwald[5] claim that fine grid relaxationis most effective after the smoothness of the increment is accounted for.Interpolation for the FMG step is bilinear. Brandt and Greenwald[5] argue that fortime-dependent problems higher order interpolation only decreases high frequency error,which tends to dissipate in parabolic systems anyway. This suggestion is utilized here,since experiments with cubic interpolation did not produce any reduction in the numberof multigrid iterations to achieve a given tolerance.Another recommendation of Brandt and Greenwald[5] is to avoid the final smoothingpass at each time step, in order to reduce aliasing from Red-Black Gauss-Seidel relaxation.Aliasing is not a major source of error in the model problem, so this suggestion does notimprove accuracy or reduce the number of fine grid iterations. For these reasons, thissuggestion was not utilized.For this problem we use a residual test with a tolerance TOL to determine the numberof fine grid iterations to perform at each time step. In the next section we show that thechoice of time-stepping scheme can affect the number of fine grid iterations to achieve aChapter 4. Convection-Diffusion Experiments^ 59given residual tolerance.4.3.3 Relative Efficiencies of Several IMEX SchemesThe model problem (4.40) was approximated using step sizes h = -121 8 and k = 0.00625and residual tolerance TOL=0.003. After the first time step, high frequency informationHF(x,y). 0.005 cos [27r(64x 63y)]was added to each of u and v, to represent the type of high frequency information thatmight be produced during a computation.For several second order IMEX schemes, the average number of fine grid iterationsat each time step was computed. The result using V-cycles is graphed in Figure 4.32.Strongly damping schemes such as SBDF and modified CNAB require roughly 1 iterationper time step. Weakly damping schemes required far more effort to solve the implicitequations accurately, because lingering high frequency components necessitate more workon the finest grid. This is evident since CNAB requires more than 2 iterations per timestep and CNLF required 3.Using a W-cycle improves the relative efficiencies of CNLF and CNAB. Even for thesecases, however, nearly twice the number of fine grid iterations were required than for morestrongly damping schemes such as SBDF and modified CNAB. See Figure 4.33 for theseresults.Even for a smaller viscosity, v = 0.02, and a much coarser mesh h = th, the per-formance of CNLF suffers. It uses about 30% more iterations than SBDF or modifiedCNAB to achieve the desired tolerance. This result uses TOL=.009 and is plotted inFigure 4.34.Thus we conclude that use of a poorly damping IMEX scheme such as CNAB orespecially CNLF can necessitate extra iterations on the finest grid in multigrid solutionsChapter 4. Convection-Diffusion Experiments^ 60to the implicit equations for small mesh Reynolds numbers, R < 2. For large meshReynolds, R > 2, this effect was not observed.4.4 Choice of IMEX SchemeBased on observations in this and previous chapters we provide a few guidelines forselecting IMEX schemes for convection-diffusion problems.4.4.1 Finite DifferencesFor the finite difference case, begin by determining an estimate for the inverse of themesh Reynolds number, h, where v represents viscosity, a and h represent characteristicspeed and grid spacing.For problems where fri << aapplication of CNLF, or a third or fourth order schemeis reasonable. A study to determine the optimal higher order scheme for ^< would61Chapter 4. Convection-Diffusion Experiments^ 62be interesting. Third and fourth order SBDF, can be applied to these problems, sincethe fl-axis is included in the stability region even though these methods were selectedprimarily for their large v properties. Of the methods we have considered in detail, CNLFhas the mildest time step restriction and is non-dissipative while third order SBDF is themost dissipative. All other second order schemes should be avoided in this case. Explicitschemes should also be considered for these problems.For 0 < Jai < -1 use of CNLF or third order SBDF appears appropriate. CNLF hasthe mildest time step restriction, but accuracy concerns could make third order SBDFcompetitive. Avoid use of SBDF in this case. For problems of this type, a study todetermine when explicit schemes are competitive would be interesting.For c''7, a theory predicts that many second order IMEX schemes have similartime step restrictions. A study to determine the method with the smallest truncationerror would be useful in this case. For greater accuracy, third order SBDF appears tobe more useful than fourth order SBDF, since its time step restriction is less severe. Ifstrong decay of high frequency spatial modes is a desirable characteristic modified CNABcan be used.For (*, > 2 , use of SBDF permits the largest stable time steps. Modified CNAB canalso be applied to problems of this type, although its time step restriction is somewhatstricter. Third order SBDF is recommended when high accuracy is needed. Numericalexperiments in Section 4.3 demonstrate that in multigrid solutions to the implicit equa-tions, application of a strongly damping method is prudent. Modified CNAB or SBDFshould be useful in such problems. Avoid use of CNLF in this case.4.4.2 Spectral MethodsBecause the eigenvalues for spectral methods are different than for finite differences as istheir meaning (see Trefethen[20]), we cannot expect the stability time step restrictionsChapter 4. Convection-Diffusion Experiments^ 63from Sections 2.7 and 3.6 to hold quantitatively. For this reason, an analysis of thelinear advection-diffusion equation for Chebyshev and Fourier spectral methods wouldbe interesting.Numerical experiments for the Burgers equation were made for small to moderatemesh Reynolds numbers. For these problems, CNLF has a very severe time step restric-tion. These computations also suggest that SBDF has the mildest time step restrictionof the methods considered. This result is particularly pronounced in the Chebyshevcollocation case.The large mesh Reynolds number case for spectral collocation was not considered.In this case, a comparison of the relative efficiencies of IMEX schemes and fully explicitschemes would be interesting.Third order SBDF appears to be an efficient method for problems where high accuracyis needed.In problems where aliasing occurs, a strongly damping scheme, such as SBDF, modi-fied CNAB or third order SBDF can be used to inexpensively reduce aliasing. Applicationof a weakly damping scheme such as CNAB or CNLF in poorly spatially resolved, aliasedcomputations should be avoided.Chapter 5Reaction-Diffusion ExperimentsThe previous chapters focus on the application of IMEX schemes to convection-diffusionproblems. We now examine properties of IMEX solutions to reaction-diffusion applica-tions.The next section introduces the prototype two chemical reaction-diffusion problem.Subsequent sections develop theory for computing this problem and demonstrate stabilityresults for a simple application.5.1 Reaction-Diffusion SystemsWe now consider the general two chemical reaction-diffusion systemu t = f (u, v) + vi Au , 1)1 > 0vt = g (u, v) + v2 0v , v2 > 0^(5.41)where the scalar functions u and v represent chemical concentrations and f and g arenonlinear. Generalization to larger systems is straightforward.To approximate the system (5.41) using an explicit algorithm frequently necessitatesunreasonably small time steps, due to the diffusive terms. An implicit algorithm, however,must contend with the expense and complexity of solving a nonlinear system at each timestep. We hope to avoid these problems by treating the nonlinear terms explicitly and thediffusive terms implicitly, using IMEX schemes.64Chapter 5. Reaction-Diffusion Experiments^ 65To demonstrate how IMEX solutions to the system (5.41) differ from those of convection-diffusion, linearize (5.41) about the steady state (u o , vo ) to obtain?it = fu(vo, vo)( 11 ) Muo, vo)(i3) viAu + 0(u 2 ) + 0(ftV) Q(v2)Vt^gu (uo , vo)(it) g, (uo, vo)( 75 ) ii20v + 0(u 2 ) + 0(itij) + 0(1) 2 )where it = u — u o and ii = v — vo . The matrix[fu(vo, vo) iv(vo, vo)gu(uo, vo) Muo, vo) or simply[fu fygu(5.42)has eigenvaluesfu gv V(fu gv) 2 4fvgu — 4fugvA1 ,2 =2Thus an IMEX approximation to (5.41) will most likely require the explicit treatment ofeigenvalues with nonzero real parts. For example, in morphogenesis applications,Re(4 2) < 0 (5.43)for the possibly complex eigenvalues, '1,2, because two of the conditions for generatingspatial patterns by a two chemical mechanism (see Murray[16]) arefu + gv < 0,fug. - fugv < 0.The spatial patterns which can form in these examples arise from the coupling ofdiffusion and reaction terms. See Murray[16] or the pioneering work by Turing[21] fordetailed derivations of this mechanism.The next section examines how to integrate equations for which the property (5.43)holds.Chapter 5. Reaction-Diffusion Experiments^ 665.2 Treatment of Reaction TermsTo determine methods for the efficient integration of reaction-diffusion problems, we con-sider the explicit case v 1 = v2 = 0. Stability properties of the explicit extrapolation areexamined for first, second and third order SBDF as well as CNLF and CNAB. Meth-ods are sought which possess large stability regions which can accommodate eigenvalueswith negative real parts. For this purpose, stability contours for the prototype equation= (i# a)x are derived as in Section 2.10. These are displayed in Figures 5.35 to 5.37.The contours for leap frog are not displayed because only values on the 0-axis arestable. The corresponding IMEX scheme, CNLF, should not be used for reaction-diffusionbecause of leap frog's inability to deal with nonzero real eigenvalues.The contours for second order Adams-Bashforth are displayed in Figure 5.35. Thestability region for this scheme includes values which are not purely imaginary. Thus thecorresponding second order IMEX schemes (e.g. CNAB and modified CNAB) are betterchoices for computing reaction-diffusion problems than CNLF.The contours for the explicit extrapolation used in SBDF are plotted in Figure 5.36.Since this stability region is significantly larger than the region for second order Adams-Bashforth, larger stable time steps should be possible when applying SBDF to reaction-diffusion problems.The contours for the explicit extrapolation used for third order SBDF are plotted inFigure 5.37. This region is smaller than for the second order case, so we expect third orderSBDF to possess a more severe time stepping restriction. This stability region containsthe 9-axis, however, so eigenvalues near the 9-axis may be more efficiently handled bythis method.Finally we consider Forward Euler, the extrapolation for all first order methods.Applying this scheme to x = (a + i#)x yields a stability region consisting of all a and #Chapter 5. Reaction-Diffusion Experiments^ 67such that11 + ii3k + aid < 1which implies that1^1( ct + _12 + #2 ___'^k i^k2 . (5.44)Thus the stability region is bounded by the circle, centre (a, /3) = (-1/k, 0), radius 1/k.Since this region is larger than any of the methods described above, we expect first orderSBDF to allow the largest stable time steps for reaction-diffusion problems.Throughout this section, stability effects from the viscous term are neglected. Toverify the usefulness of this simplified stability analysis, numerical experiments for areaction-diffusion problem are carried out in the next section.Chapter 5. Reaction-Diffusion Experiments^ 68Chapter 5. Reaction-Diffusion Experiments^ 695.3 Computations5.3.1 The Schnackenberg SystemFor our model problem we consider the one-dimensional Schnackenberg system (seeSchnackenberg[18]),ut = d(a — u u 2v) uxs^vt = d(b — u2v) vvss^ (5.45)having periodic boundary conditions on the interval [0,1].This model has been used by Murray[16] in the context of morphogenesis. To considerpattern formation consistent with morphogenesis we use the initial conditionscos(2irjx)u(x, 0) = a + b 0.0004J=1)v(x, 0)^cos (27rjx^ + 0.0004 E^(a + b)2^j=1These initial perturbations from the steady state(U,^= (a b, (a + b)2)represent the random modal frequencies that could be present in a biological system.We select parametersa = 0.2b = 2.0andd = 300.Chapter 5. Reaction-Diffusion Experiments^ 70From Murray[15], spatial patterns can be generated using these parameter values. Havingselected a, b and d, an upper bound for the time step size assuming no diffusion can bedetermined. The Jacobian matrix (5.42) evaluated at (u, has approximate eigenvalues—603 ± 268i for these parameters. Using (5.44), the time step size, k, will have to betaken less than about 0.0027 to ensure that the eigenvalues lie inside the stability regionfor first order SBDF at t = 0.Computations are carried out for 60 < v < 160. For v < 60, spatial error overwhelmstime stepping error because patterns form very slowly or not at all. For v > 160, thequalitative behaviour of the final patterns is similar to that of v = 160.There are two qualitatively different solutions. In Figures 5.38 and 5.39, results forv = 60 are provided for u and v. In Figures 5.40 and 5.41, the results for v = 160 aredisplayed. These computations are accurate to within 1% of the exact solution and usefirst order SBDF with fourth order centred differences spatially. The discretization stepsizes are k = 0.0005 in time and h = 0.05 in space.In the next section, several stability properties for IMEX schemes are demonstratedfor the system (5.45).5.3.2 One-Dimensional ComputationsTo illustrate stability properties for IMEX schemes, solutions for the model problem (5.45)were computed using discretization step sizes h = 0.05 in space and k = 0.00025 in time.Using several IMEX schemes, computations are performed to time t = 0.8 using fourthorder centred (5-point) differences in space.The max norm relative errors for these schemes are computed by comparing results tothose for SBDF using h = 0.0125 and k = 0.00003125. The resultant errors are plottedagainst v for 100 < v < 160 in Figure 5.42. For 60 < v < 100 the errors for all methodsexcept CNLF overlap since spatial error dominates.Chapter 5. Reaction-Diffusion Experiments^ 711.4 L02.8-2.6-2.4-2.23.0.0.20.0.300*I.I.I.I.I.I,.4^.5^.6^.7^.8^.9oz.co0. 0 0.0^ft,0.iFigure 5.38: Concentration, u, at t = 0.8 for v = 606.05.55.04.54.03.53.02.59*co0.0 0.02.01.51.00.8.2^.3 .4^.5 .6 .7I^,^I .8 .9 1 0xChapter 5. Reaction-Diffusion Experiments^ 72Figure 5.40: Concentration, u, at t = 0.8 for v = 160Chapter 5. Reaction-Diffusion Experiments^ 73For 60 < v < 160, CNLF is unstable as predicted by the theory. This methodis unstable even for k = 2.5 x 10 -6 . CNAB, modified CNAB and third order SBDFproduce large errors over much of the interval. First order SBDF and SBDF appearstable over the entire interval. The error for first order SBDF is not plotted because itnearly overlaps that for SBDF.To determine the nature of the errors, computations for v when v = 160 are plotted inFigure 5.43. From Figure 5.42, the SBDF computation is nearly exact. CNAB, modifiedCNAB and third order SBDF all have incorrect amplitudes. Furthermore, the dominantmodal frequency for third order SBDF and CNAB is not that of the exact solution.Similar effects occur even for initial perturbations 100 times larger.From Figure 5.44, even for a smaller value of h, h = 1/80, the results are qualitativelysimilar. The error for modified CNAB is not plotted because it nearly overlaps that ofthird order SBDF. In practice, the value of h would not usually so small for these valuesof v because the coarser mesh, h = 1/20, yields an error of less than 1% using first orderSBDF.An increase in the time step size to k = 0.0005 leaves the error for first order SBDFessentially unchanged, whereas SBDF has an error greater than 35% for v > 120. Aspredicted by the theory, the method with the largest explicit stability region allows thelargest stable time steps.The results of this section are consistent with the hypothesis that numerical methodspossessing large stability regions in the absence of diffusion are preferred for reactiondiffusion problems. This is because even for small values of h it is possible for reactionterms to dominate diffusion if parts of the solution are linear in x. To illustrate thisbehaviour, consider the Schnackenberg system (5.45) with a = b = 0 and d = 300, subjectto periodic boundary conditions and initial conditions u(x, 0) = u0 and v(x, 0) = 0.Chapter 5. Reaction-Diffusion Experiments^ 74111-1-111111[„11111„11110 -2LOO•_SBDF10 -3Sapp.100^105^110^115^120^125^130^135^140^145^150^155^160ViscosityFigure 5.42: Relative Errors for Several IMEX Schemes for Various v and h = 20Applying second order centred differences and first order SBDF to this system yields,n+1^rt^n+1^ +1^n+1Zia ^u u •^—^•^-I- u3+ 11 "3 = 300u +rat^3-h2By symmetry, url = uni +1 = unl thus we obtain0+ 1 = (1 — 300k)03^ 3which implies the stability time step restriction of k = 150 for any h. This is the samerestriction that occurs without diffusion.The next section considers the two-dimensional case.5.3.3 Two-Dimensional ComputationsFor our next problem, we consider the two-dimensional Schnackenberg systemu t = 300(0.2 — u u 2v) + Auv t = 300(2.0 — u 2v) 1200v. (5.46)10 -2CO01^I^I10 -3 rL^0Or2 0CDU7oL10-5^SBDFSBDF SBDFChapter 5. Reaction-Diffusion Experiments^ 75.55.50.45.40.35.30.25.20.15.10.0500^ .3^.4^.5^.6^.7^.8^.9^1 0XFigure 5.43: IMEX Computations for v, at t = 0.8 for v = 160.210-6 ^la^1^i1^1,1^I^I^1^1,1^I^1^1^I^I 'tilt^1^1^I^1^I^I^1^1^I^1^1^1,1^I I^I^1^1 100^105^110^115^120^125^130^135^140^145^150^155^160ViscosityFigure 5.44: Relative Errors for Several IMEX Schemes for Various v and h =Chapter 5. Reaction-Diffusion Experiments^ 76subject to periodic boundary conditions on the interval 12 = [0,1] x [0, 1] and initialconditionsu(x, y, 0) =cos(271-j(x^y))2.2^0.001+^2_,v(x, y, 0) = 0.41322^0.001^ ov--,8^c s(21-j(x+y)) .j=.1The problem (5.46) was computed to time t = 0.2 using discretization step sizesh = 1/32 in space and k = 0.001 in time. The standard 5-point Laplacian was usedto approximate the diffusion terms. A multigrid solution to the implicit equations wasperformed using V(1,1) cycles with one iteration per time step. The components of thisiterative algorithm are described in Section 4.3.Using first order SBDF, this computation has a max norm relative error of less than0.44% when compared to a computation with SBDF using k = .0001 and h = —1218• Theplot of u for this computation is displayed in Figure 5.45. Similar to the one-dimensionalcase, SBDF, third order SBDF, CNAB and modified CNAB all produce relative maxnorm errors greater than 1 for this value of k. Refining the time step to k = 0.0005produces errors of less than 0.44% for SBDF and first order SBDF. Modified CNAB,CNAB and CNLF still produce relative errors greater than 1. All the methods exceptCNAB and CNLF have errors of less than .44% when the time step k = 0.00025 is used.The results for the refined mesh, h = 1/64, are qualitatively similar.Using the results from this and previous sections, a description of how to chooseIMEX schemes for reaction-diffusion problems is presented next.5.4 Choice of IMEX SchemeBased on the observations in this and previous chapters we provide a few guidelines forselecting IMEX schemes for reaction-diffusion problems. These ideas can also be appliedChapter 5. Reaction-Diffusion Experiments^ 77Reaction Diffusion"47:40:4:4;lb\\7"Ar!♦ val .5 Or.^4■''''..44P-1•41111r4rArArAwavAl701W■^44"4014;;411:4SAdrar:1010/41.\.\\VICAVA rA "Aft, 46.4ratitairgsfol\•1Stly41•4■4741•Ardrifrns"^ 004.41004•411PAM,-4: 414P6,4412,41IIIArlif 4.4,04nr.dur41.P.Nprars'.1,40.4,4r AiirArAgg,,A41,48.4,,,Wargir^ 4IP04.3.1. 4•47'•'-argar" 101Figure 5.45: Concentration, u, at t = 0.2 for v = 120Chapter 5. Reaction-Diffusion Experiments^ 78to problems where explicitly treated terms have eigenvalues with negative real partswhose size is not small compared to the implicitly treated terms.The primary observation is that CNLF should not be used in reaction-diffusion prob-lems. This scheme cannot accommodate the real eigenvalues of the reaction term, andwill likely become unstable.SBDF methods are excellent choices for reaction-diffusion problems. These methodsallow large time steps in comparison to other methods of the same order because theirexplicit extrapolations include more of the negative a-axis than other schemes of thesame order.Finally, application of an IMEX scheme of lower order than the spatial discretizationshould be considered. Lower order methods allow larger time steps because their explicitextrapolations tend to possess larger stability regions.Chapter 6ConclusionsFrom the analysis and experiments of the previous four chapters, several conclusions andrecommendations for future work are made. These are given in the next two sections.6.1 ResultsThe analysis of the general, linear multistep IMEX scheme demonstrates that s-stepIMEX schemes can have at most order s accuracy, and that such schemes form an s-parameter family of methods. This fact leads us to consider only s-step, order s schemesthroughout this thesis.A parameterization for first order, one-step schemes was given. Of these methods,first order SBDF (2.15) was recommended because it allows large stable time steps forv > 1/k and produces a strong asymptotic decay of high frequency spatial modes.A parameterization for second order, two-step schemes was given. Included in thistwo-parameter family of schemes are CNLF (2.19), CNAB (2.17), SBDF (2.20) and ournewly proposed scheme, modified CNAB (2.18). An analysis for the linear advection-diffusion equation showed that for small mesh Reynolds numbers (2.24), R < 2, SBDFproduces the strongest asymptotic decay of high frequency spatial modes and allows thelargest stable time steps. CNLF was shown to allow the largest stable time steps forlarge mesh Reynolds numbers, R > 2. The modified CNAB scheme is recommended forthe small mesh Reynolds number case when a small truncation error along with a strongdecay of high frequency error modes is desired.79Chapter 6. Conclusions^ 80For the third order case, a three-parameter family of three-step schemes was given.It was shown that third order SBDF (3.26) produces the strongest asymptotic decay ofthese schemes and that it is stable for all v > 0. For any v > 0, however, one of CNLF orSBDF allows larger stable time steps. A fourth order scheme, fourth order SBDF (3.36)was also given. This scheme was shown to be stable for all v > 0. It possesses an evenmore severe time-stepping restriction than third order SBDF.Computations for solutions to a variable coefficient problem were carried out usingsecond order, centred differences spatially. Several of the previously identified theoreticalresults were demonstrated. In particular, SBDF allowed the largest stable time stepsfor small mesh Reynolds numbers and CNLF allowed the largest stable steps for largemesh Reynolds numbers. The errors for CNAB and modified CNAB were very similar,supporting the claim that these methods have similar truncation errors. Stability formodified CNAB was also similar to CNAB. Third order SBDF was shown to possessa more severe time-stepping restriction than many second order methods. Nonetheless,this method is recommended when high accuracy is needed because it produces a muchsmaller error than second order schemes when stable.Based on solutions for a nonlinear problem we conclude that CNLF, third and fourthorder SBDF perform well for the zero viscosity case. Of these schemes, CNLF was non-dissipative while third order SBDF was the most strongly dissipative. Other second orderschemes should be avoided in this case.Using Chebyshev and Fourier spectral collocation, solutions to the Burgers equa-tion for small to moderate mesh Reynolds numbers were computed using various IMEXschemes. For these problems, CNLF had a very severe time step restriction. These com-putations also suggest that SBDF had the mildest time step restriction of the methodsconsidered. Although it had a more severe time-stepping restriction than many secondorder schemes, third order SBDF appeared to be useful when high accuracy was required.Chapter 6. Conclusions^ 81Further pseudospectral computations for the Burgers equation demonstrated that astrongly damping time-stepping scheme such as SBDF, modified CNAB or third orderSBDF can be used to inexpensively reduce aliasing. Application of a weakly dampingscheme such as CNLF or CNAB in poorly spatially resolved, aliased computations isdiscouraged.Using a multigrid method to solve the implicit equations, solutions to the 2D convection-diffusion equation for small mesh Reynolds numbers were computed using various IMEXschemes. These results clearly demonstrate that use of a strongly damping time-steppingscheme may lead to fewer fine grid iterations in multigrid than use of a weakly damp-ing scheme. In particular, use of SBDF or modified CNAB produced a more efficientmultigrid method than did use of CNLF or CNAB. The idea of using a strongly dampingtime-stepping scheme to reduce the number of iterations in multigrid should also extendto fully implicit schemes.The reaction-diffusion problem was considered. For this problem, it is demonstratedthat schemes which possess large stability regions in the absence of diffusion are preferred.Such schemes are best able to cope with the growth associated with the reaction terms.Of the IMEX schemes considered, first order SBDF allowed the largest time steps andsecond order SBDF allowed the largest stable time steps of second order schemes. Useof CNLF is strongly discouraged in reaction-diffusion problems.Perhaps the most surprising conclusion of this study is the poor showing of the mostpopular scheme currently in practice, namely CNAB. In particular, our newly proposedmodified CNAB scheme (2.18) appears to supersede CNAB because it almost alwaysperformed as well as, or much better than, CNAB.Chapter 6. Conclusions^ 826.2 Future WorkThere are many opportunities for future work in IMEX schemes.An investigation to determine when IMEX schemes are more efficient than fully ex-plicit or fully implicit schemes would be of interest. For the convection-diffusion case ananalysis to determine when a fully explicit treatment is competitive is especially needed.For this case, we expect fully explicit schemes to be more efficient for very large meshReynolds numbers because the time-stepping restriction arises primarily from the con-vection term.A study of the constant of the truncation error for IMEX schemes would be veryinteresting, especially for the second order case. Such an analysis would help determinewhich second order scheme to apply to a particular problem. However, this is far lesssimple than for usual ODE multistep schemes.It is clear that algorithms using IMEX schemes are easier to implement than thoseusing fully implicit schemes in reaction-diffusion problems. We have not considered howefficient IMEX schemes are in comparison to fully implicit schemes for these systems. Astudy making this comparison would be of interest.The analysis of this thesis has focussed on the finite difference case. Because theeigenvalues for spectral methods are different than for finite differences as is their mean-ing, we suspect that parts of the theory will not apply to spectral methods. For thisreason, an analysis of the linear advection-diffusion equation for Chebyshev and Fourierspectral methods would be useful.Future work in time-dependent multigrid is also needed. This thesis has only ex-amined the small mesh Reynolds number case. The large mesh Reynolds number caseshould be examined because it may also yield interesting behaviour.Throughout this work, periodic boundary conditions were assumed. In practice, weChapter 6. Conclusions^ 83usually have to impose non-periodic boundary conditions. An analysis to determine sta-ble methods of imposing Dirichlet and Neumann boundary conditions for IMEX schemeswould be very useful.This thesis has applied a one-step scheme, first order SBDF, with a very small timestep to obtain the starting values in multistep IMEX schemes. A study to determinesimple, efficient methods for obtaining these initial values would be of interest.Further analysis of third and fourth order schemes is also possible. For example, itwould be desirable to determine if schemes exist which allow significantly larger stabletime steps than third or fourth order SBDF schemes.For systems which are strongly convective, the phase lag can dominate the error. Ananalysis of this source of error for IMEX schemes should also be carried out.Finally, we have only considered the constant time step size case. In practice, variabletime-stepping algorithms are often preferred. A study to examine how to implementvariable time-stepping for various IMEX schemes would be interesting.Bibliography[1] C. Basdevant et al. Spectral and finite difference solutions of the Burgers equation.Computational Fluids, 14:23-41, 1986.[2] J.P. Boyd. Chebyshev 1 Fourier Spectral Methods. Springer-Verlag, 1989.[3] M.E. Brachet et al. Small-scale structure of the Taylor-Green vortex. J. FluidMechanics, 130:411-452, 1983.[4] A. Brandt. Guide to multigrid development. In W. Hackbusch and U. Trottenberg,editors, Multigrid Methods, pages 220-312. Springer-Verlag, Berlin, 1982.[5] A. Brandt and J. Greenwald. Parabolic multigrid revisited. In Proc. 3rd EuropeanConference on Multigrid Methods, Bonn, Oct. 1990.[6] C. Canuto, M.Y. Hussaini, A. Quarteroni, and T.A. Zang. Spectral Methods in FluidDynamics. Springer-Verlag, 1987.[7] Bruce W. Char et al. MapleV Language Reference Manual. Springer-Verlag, 1991.[8] C.W. Gear. The automatic integration of stiff ordinary differential equations. InProc. IFIP Congress, pages A81-85, Ljubjana, Yugoslavia, 1986.[9] G.H. Golub and C.F. Van Loan. Matrix Computations. Johns Hopkins Univ. Press,Baltimore, second edition, 1989.[10] W. Hackbusch. Multi-Grid methods and Applications. Springer-Verlag, 1985.[11] E. Hairer, S.P. Norsett, and G. Wanner. Solving Ordinary Differential Equations I.Springer-Verlag, 1987.[12] C. Hirsch. Numerical Computation of Internal and External Flows. Wiley, 1988.[13] J. Kim and P. Moin. Application of a fractional-step method to incompressibleNavier-Stokes equations. J. Computational Phys., 59:308-323, 1985.[14] J.J.H Miller. On the location of zeros of certain classes of polynomials with applica-tions to numerical analysis. J. of the Institute of Mathematics and Its Applications,8:397-406, 1971.84Bibliography^ 85[15] J.D. Murray. Parameter space for Turing instability in reaction diffusion mecha-nisms: a comparison of models. J. Theoretical Biology, 98:143-162, 1982.[16] J.D. Murray. Mathematical Biology. Springer-Verlag, 1989.[17] R. Peyret and T. Taylor. Computational Methods for Fluid Flow. Springer-Verlag,1983.[18] J. Schnackenberg. Simple chemical reaction systems with limit cycle behaviour. J.Theoretical Biology, 81:389-400, 1979.[19] J.C. Strikwerda. Finite Difference Schemes and Partial Differential Equations.Wadsworth & Brooks/Cole, 1989.[20] 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, pages237-253. Clarendon Press, 1988.[21] A.M. Turing. The chemical basis of morphogenesis. Phil. Trans. Roy. Soc. Lond.,B237:37-72, 1952.[22] J.M. Varah. Stability restrictions on second order, three level finite differenceschemes for parabolic equations. Siam J. Numerical Analysis, 17(2):300-309, 1980.[23] R. Varga. Matrix Iterative Analysis. Prentice-Hall, Englewood Cliffs, NJ, 1962.
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Implicit-explicit methods for time-dependent PDE’s
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Implicit-explicit methods for time-dependent PDE’s Ruuth, Steven J 1993
pdf
Page Metadata
Item Metadata
Title | Implicit-explicit methods for time-dependent PDE’s |
Creator |
Ruuth, Steven J |
Date Issued | 1993 |
Description | Various methods have been proposed to integrate dynamical systems arising from spatially discretized time-dependent partial differential equations. For problems with terms of different types, implicit-explicit (IMEX) schemes have been used, especially in conjunction with spectral methods. For convection-diffusion problems, for example, one would use an explicit scheme for the convection term and an implicit scheme for thediffusion term. Reaction-diffusion 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 aliasing reduction for spectral methods. For the prototype linear advection-diffusion equation, a stability analysis for first, 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 identified. Numerical experiments demonstrate that weak decay of high frequency modes can lead to extra iterations on the finest grid when using multigrid computations with finite difference 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 findings are demonstrated on several examples. |
Extent | 3517774 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
FileFormat | application/pdf |
Language | eng |
Date Available | 2008-09-15 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0079818 |
URI | http://hdl.handle.net/2429/1950 |
Degree |
Master of Science - MSc |
Program |
Mathematics |
Affiliation |
Science, Faculty of Mathematics, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 1993-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-ubc_1993_spring_ruuth_steven.pdf [ 3.35MB ]
- Metadata
- JSON: 831-1.0079818.json
- JSON-LD: 831-1.0079818-ld.json
- RDF/XML (Pretty): 831-1.0079818-rdf.xml
- RDF/JSON: 831-1.0079818-rdf.json
- Turtle: 831-1.0079818-turtle.txt
- N-Triples: 831-1.0079818-rdf-ntriples.txt
- Original Record: 831-1.0079818-source.json
- Full Text
- 831-1.0079818-fulltext.txt
- Citation
- 831-1.0079818.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0079818/manifest