Spatially dispersionless, unconditionally stable FCAD solvers for variablecoefficient PDEs
Abstract
We present fast, spatially dispersionless and unconditionally stable highorder solvers for Partial Differential Equations (PDEs) with variable coefficients in general smooth domains. Our solvers, which are based on (i) A certain “Fourier continuation” (FC) method for the resolution of the Gibbs phenomenon, together with (ii) A new, preconditioned, FCbased solver for twopoint boundary value problems (BVP) for variablecoefficient Ordinary Differential Equations, and (iii) An Alternating Direction strategy, generalize significantly a class of FCbased solvers introduced recently for constantcoefficient PDEs. The present algorithms, which are applicable, with highorder accuracy, to variablecoefficient elliptic, parabolic and hyperbolic PDEs in general domains with smooth boundaries, are unconditionally stable, do not suffer from spatial numerical dispersion, and they run at FFT speeds. The accuracy, efficiency and overall capabilities of our methods are demonstrated by means of applications to challenging problems of diffusion and wave propagation in heterogeneous media.
Keywords: Highorder methods, Alternating Direction Implicit schemes, numerical dispersion, variable coefficient problems.
1 Introduction
We present fast, spatially dispersionless and unconditionally stable highorder solvers for Partial Differential Equations (PDEs) with variable coefficients in general smooth domains. Our algorithms, which generalize significantly a class of solvers [bruno10, lyon10] introduced recently for constantcoefficient PDEs, are based on (i) A certain “Fourier continuation” (FC) method [bruno10] for the resolution of the Gibbs phenomenon, together with (ii) A new, preconditioned, FCbased solver for twopoint boundary value problems (BVP) for variablecoefficient Ordinary Differential Equations (ODE), and (iii) The Alternating Direction Implicit (ADI) methodology [peaceman55, douglas56]. One of the main enabling elements in our overall FCAD algorithm (FourierContinuation AlternatingDirections) is the new solver (ii) for twopoint boundary value problems with variable coefficients. Relying on preconditioners that result from inversion of oversampled finitedifference matrices together with a new methodology for enforcement of boundary conditions and the iterative linear algebra solver GMRES, this algorithm produces rapidly the solutions required for FCAD timestepping in the variablecoefficient context. (A nonoversampled finitedifference preconditioner related to but different from the one used here was introduced in [Sun96] in the context of orthogonal collocation methods for equations with constant coefficients.) The resulting PDE solvers, which in practice are found to be unconditionally stable, do not suffer from spatial numerical dispersion and they run at a computational cost that grows as with the size of the computational grid. A variety of examples presented in this paper demonstrate the accuracy, speed and overall capabilities of the proposed methodology.
The variablecoefficient FCAD algorithms introduced in this paper enjoy all the good qualities associated with the constant coefficient solvers presented in [bruno10, lyon10]: the performance of the new solvers compare favourably, in terms of accuracy and speed, with those associated with previous approaches; a detailed discussion in these regards can be found in the introductory sections of [bruno10, lyon10]. In particular, in this paper we demonstrate the highorder, essentially dispersionless character of the new FCAD solvers by means of solutions to parabolic and hyperbolic problems. For example, the results presented in Section 7.3.1, which include FCAD fixedaccuracy solutions at fixed numbers of pointsperwavelength for problems of sizes ranging from one to onehundred wavelengths in size, demonstrate the spatial dispersionlessness of the FCAD algorithm for hyperbolic equations. Further, as shown in Table 1, for example, the proposed FCAD algorithm can evolve a solution characterized by onemillion spatial unknowns in a computing time of approximately seconds per time step in a singlecore run.
The remainder of this paper is organized as follows: after the introduction in Section 2.1 of the variablecoefficient PDEs we consider, Section 2.2 details the ADI approximations we employ. Section 3 presents the FC method, including, in Section 3.3, a special version of the FC algorithm we need to tackle variablecoefficient differential equations. Sections 4 and 5 then describe our new FCbased solver for twopoint BVP, which include 1) Iterative FCbased solvers for periodic ordinary differential equations in a certain “continued” periodic context, and 2) Techniques that enable enforcement of boundary conditions in presence of either sharp or diffuse boundary layers. The combined FC/ADI scheme is described in Section 6. The overall properties of the resulting FCAD PDE solver, finally, are demonstrated in Section 7 through a variety of numerical results.
2 Preliminaries
2.1 Variable coefficients PDEs
This paper presents FCbased solvers for linear equations containing timeindependent but spatially variable coefficients. While the methods we present are applicable to any partial differential equation for which an ADI splitting is available, for definiteness we focus on the basic variablecoefficient parabolic and hyperbolic problems
(1) 
and
(2) 
in a bounded open set with smooth boundary and within the time interval . Here , , , , and are suficiently regular functions defined in ; additionally the coefficient functions and are assumed to satisfy the coercivity conditions
for some positive constants , , and , while the initial and source functions are required to verify the relevant compatibility conditions in , namely, for the parabolic problem (1) and , for the hyperbolic problem (2).
2.2 Alternating Direction schemes
This section presents the ADI splitting schemes we use for the solution of the PDE problems (1) and (2). These splitting schemes result as adequate generalizations of the PeacemanRachford scheme [peaceman55] for diffusion problems and the noncentered scheme [lyon10] for the wave equation to the present variablecoefficient context. In what follows denotes the time step used in the computational time interval ; it is assumed that for a certain positive integer . Letting for integer and even fractional values of (e.g., ), we have, in particular, .
2.2.1 Diffusion equation
Given the initial values , the right handsides and , and the boundary values , the exact solution for of the diffusion problem (1) satisfies the CrankNicolson relation [peaceman55, bruno10]
(3) 
Following [peaceman55], an ADI scheme can be obtained from the CrankNicolson iteration: denoting by the corresponding approximation of (for integer values ) and using the intermediate quantity , the ADI scheme is embodied in the equations
(4) 
(with boundary condition on ) and
(5) 
(with boundary condition on ); cf. [bruno10].
An algorithm based on the ADI iteration (4)(5) can be conveniently obtained as a sequence of four operations involving two additional auxiliary quantities and :

Initialize and as
Then, for ,

Obtain by solving the boundary value problem
(6) 
Update according to
and, finally,

Obtain by solving the boundary value problem
(7)
The ADI scheme thus requires solution of the onedimensional boundary value value problems (D2) and (D4) (see Section 4) and solution updates (D1) and (D3). Each one of these operations involves differential operators with respect to a single spatial variable—as it behooves an ADI discretization [marchuk90].
Remark 1.
An implementation of this algorithm for the case in which the coefficients and are constant was put forth in [bruno10]. In that reference it was noted that, for , is a globally second order accurate approximation of the exact solution of the diffusion equation (1): the error at any fixed time step is a quantity of order . This is in spite of the approximation , which is necessary to enable applicability to complex domains, and which induces a secondorder local truncation error in step (D2); see [bruno10, Remark 3.2] for details. Numerical experiments we present in this paper reveal once again a secondorder global error in the solutions resulting from the scheme above. As shown in reference [bruno10] and Section 7.2, further, solutions of higher order of temporal accuracy can be extracted from the solutions produced by the ADI scheme above by means of the Richardson extrapolation methodology.
2.2.2 Wave equation
To derive our Alternating Direction scheme for the linear wave problem (2), in turn, we follow [lyon10] and note that the exact solution satisfies the discrete relation
(8) 
where, once again, and . Following [lyon10] we split the stiffness term and thus obtain the ADI timestepping scheme

Initialize and as
Then, for ,

Obtain by solving the boundary value problem,
(9) 
Obtain as the solution of the boundary value problem
Remark 2.
Note that the expressions (D2)(D4), and (W2)(W3) of the alternatingdirection ODEs for the heat and wave equations incorporate a division by the lowestorder variablecoefficient in equations (3) and (8). This is an essential detail of our algorithm: if such a division by is not incorporated in the algorithm, the iterative GMRES solution of these ODEs (which is presented in Section 5.1) would require large numbers of iterations. In the divided form, in contrast, the matrix of the linearalgebra problem is close to the identity for small , and small numbers of GMRES iterations suffice to yield highly accurate ODE solutions. In fact, we have found in practice that the divided ODE forms give rise to small numbers of iterations even for large values of .
Remark 3.
It is easy to check that the scheme (W1)(W3) is first order consistent. We refer to [lyon10, Section 5] for a discussion of the global order of accuracy of the algorithm; in practice, and in agreement with that reference, we find the algorithm produces solutions with global first order accuracy. As shown in reference [bruno10] and Section 7.3.2, further, solutions of higher order of temporal accuracy can be extracted from the solutions produced by the ADI scheme above by means of the Richardson extrapolation methodology.
Problems (D2), (D4), (W2) and (W3) amount to twopoint BVP of the form
(10)  
(11) 
where the righthand side and the variable coefficients and are bounded smooth functions defined in the interval , with for some constant . More precisely, the ODE coefficients and the righthand sides of the problems (D2) and (D4), (W2) and (W3) are given by
(12) 
where for (D2) and (D4)
(13) 
while for (W2) and (W3)
(14) 
To solve the twopoint BVP (10)(11) numerically we utilize a discrete method based on three main elements: 1) The Fourier Continuation method (see Section 3) to produce accurate Fourierseries approximations of the ODE source terms and variable coefficients; 2) A specialized Fourier collocation method for solution of ODE boundary value problems (see Section 4.1)—which, in view of item 1), can be applied to nonperiodic boundaryvalue problems without the accuracy degradation associated with the Gibbs phenomenon; and 3) A new strategy for the enforcement the boundary conditions in the context arising from items 1) and 2) (see Sections 4.2 and 4.3).
3 Fourier Continuation and scaled FC(Gram)
3.1 Discrete Fourier analysis background
Let denote the space of times continuously differentiable periodic functions of period . The discrete Fourier series of is given by
Here for even and for odd,
(15) 
denotes the dimensional space of trigonometric polynomials, and the amplitudes are given by the trapezoidalrule expression
Note that, for conciseness, the “degree” is not explicitly displayed in the notation .
We point out that, as is wellknown [hesthaven07, boyd01],

An element of the set is determined uniquely by its values at the equispaced grid , where
(16) 
The discrete Fourier operator is an interpolation operator from into , that is, for all .
3.2 Fourier Continuation and FC(Gram)
The Fourier Continuation method [bruno10, albin11] is an algorithm which, acting on a set of values of a function at equidistant points , produces a trigonometric polynomial , of a given degree and on a given interval with , which approximates the function closely on the interval . Note that the endpoints and are not required to belong to the Fourier grid . Clearly, for , the trigonometric interpolation operator leads to a naive Fourier Continuation procedure which gives rise to the well known Gibbs ringing near and —unless is a smooth periodic function of period . The selection of a larger expansion interval allows for a smooth transition between the values near to the values near , and thus enables highly accurate Fourier approximation in the interval , avoiding the undesirable Gibbs ringing effect and related accuracy deterioration.
A number of algorithms has been introduced which provide accurate Fourier continuations for smooth nonperiodic functions (including, for example, the FC(SVD) method [bruno10], which relies on Singular Value Decompositions and the related algorithms [bruno07, boyd02, bruno03]). Use of the spectrally accurate Fourier continuations as a component of efficient PDE solvers, however, must rely on correspondingly efficient Fourier Continuation algorithms. For that purpose, an accelerated FC method, the FC(Gram) procedure, was introduced recently [bruno10, lyon10, albin11]. In this section, we present a brief description of FC(Gram) method; full details can be found in [bruno10]. A new “scaled” version of the FC(Gram) approach, which is necessary for the applications presented in this paper, is introduced in Section 3.3.
For the present description of the FC(Gram) method, without loss of generality we assume , as in references [bruno10, lyon10]. The FC(Gram) algorithm is based on the leftmost and rightmost grid values of within the subintervals and , with . Fixed the number of extension points (belonging to with ) where has to be computed, the continuation problem is reduced to seek a smooth periodic function in . For adequate reference note that our quantities , and are denoted by , and in the contribution [bruno10].
To compute , the leftmost grid values of in are mapped in . Then, they are projected on certain Gram bases of orthogonal polynomials—with orthogonality dictated by the natural discrete scalar product defined by the grid points. Finally, using the orthogonal polynomial expansion, a periodic matching function , in the interval is computed. This function consists of a sum of contributions of each polynomial projection (up to degree ) and blends smoothly the values at the left and right subintervals, and . The grid values of in give the required grid values of . This procedure can be efficiently implemented by precomputing the sets of matching functions and associated to the even and odd pairs of the polynomial basis in the interval and then, each FC(Gram) continuation in an arbitrary interval can be obtained by using an affine transformation (see [bruno10, Section 2.3] for details).
It is clear that the accuracy of the Fourier continuation operator depends on the number of Fourier grid points contained in , the maximum degree of polynomials used in the polynomial projections, the number of right and left grid points, and also on the number of extension points outside the interval (see item 2 in Section 7 for an indication on the actual values of these parameters used in our numerical examples).
3.3 Scaled FC(Gram) algorithm
The FCbased PDE solvers presented in this paper rely on the FC methodology as a tool for solving twopoint BVPs for ordinary differential equations with variable coefficients. As discussed in Section 5.3, to obtain an accurate solution of the relevant twopoint BVPs for a given PDE problem in the present variablecoefficient context, our method requires, unlike previous FCbased approaches, use of numbers of extension points and associated extension intervals of length that vary (linearly) with (see Section 5.2). Use of such dependent extension intervals in the procedure described in Section 3.2 entails evaluation of correspondingly dependent sets of matching functions and —a procedure that is inelegant and computationally expensive.
To overcome this difficulty, a slightly modified “scaled” version of the FC(Gram) algorithm is introduced in this section. This approach is based on use of a set of precomputed matching functions , as described in Section 3.2, in a fixed interval and for a certain fixed value . To obtain inexpensively a matching function at a new number of extension points in a new interval (where it is assumed that ), the “scaled” procedure uses as a matching function a composition of the form
(17) 
where is a smooth diffeomorphism. To preserve the point values of in and , the scaling function is assumed to be map linearly the intervals and onto and , respectively. Throughout this paper we use the scaling function
(18) 
where, letting denote the largest integer less than or equal to , we have set , and the auxiliary function is given by
In what follows, the scaled FC(Gram) continuation of a function resulting from this procedure will be denoted by either of the following symbols
(19) 
The first two of these notations can of course be used only when the value of is either inconsequential or implicit in the context. The scaled FC(Gram) algorithm produces Fourier continuations for varying values of and at a computational cost that is essentially the same as that required by the original FC(Gram) procedure for a fixed number of extension points.
Remark 4.
Once the even and odd continuation of the Gram polynomials are computed and stored in memory, the FC(Gram) and scaled FC(Gram) continuations of any function only involve the computation of inner products of size , which requires sums and products, and the evaluation of a trigonometric interpolant, which is performed by means of an FFT of size . Since and are independent of , both the scaled and unscaled versions of the FC(Gram) algorithm run at a computational cost of operations.
4 FC BVP solver I: particular solution and boundary conditions
This section presents an FCbased method for solution of variablecoefficient BVPs of the form (10)–(11). As explained in what follows, this approach relies on the scaled FC(Gram) method to produce a periodic extension of the nonperiodic problem (10)–(11) to an interval , with boundary values at the endpoints of the original interval . As detailed in Section 4.1, a particular solution for equation (10) can easily be produced on the basis of the FC(Gram) method. The enforcement of the boundary conditions, which requires some consideration, is presented in Sections 4.2 and 4.3. In the first one of these sections a direct numerical approach is presented for the evaluation of solutions of the homogeneous problem (10)–(11) () with nonzero boundary values, which can be used to correct the boundary values of a particular solution. In Section 4.3 a complementary approach is introduced, which can effectively treat challenging boundary layers and stiff equations that arise as small time steps and correspondingly small coefficients are used (cf. equations (6), (7) and (9)).
4.1 FCbased particular solution
To obtain a periodic embedding to the interval of the BVP (10)–(11) we utilize Fourier continuations of the functions , and . To preserve the ellipticity of the problem, however, we must ensure that the extension used for the coefficient takes on strictly positive values. To produce such strictly positive extension of a positive function we utilize a infinitely differentiable diffeomorphism such as, for instance, . Using such a diffeomorphism we define certain “limited extensions” of a positive function by means of the expression
(20) 
where denotes the Fourier continuation procedure defined in (19). The expression (20) ensures that the values of the periodic extensions are bounded by above and below by the positive values and , respectively. Hence, the “periodic embedding” of the ODE problem (10)(11) is given by
(21) 
where, in accordance with equation (19), and denote the scaled FC(Gram) continuation of the functions and .
To produce approximate periodic solutions of period of the ODE problem (10), our algorithms relies on Fourier collocation: a numerical solution of the form
is sought that satisfies (10) at each one of the grid points , , that is
(22) 
(if problem (10)(11) has a unique solution, this system of equations is not singular, see e.g. [boyd01].) Since the solution is determined uniquely by its values for the vector of point values of , are used as the unknowns of the problem. This linear system of equations is solved by means of the iterative solver GMRES; see Section 5.1 for details.
Remark 5.
Clearly, the particular solution described in this section does not generally satisfy the boundary conditions imposed in (11). The necessary corrections are described in the following section.
4.2 Boundary conditions I: exterior sources
To enforce the necessary boundary conditions in (11) we consider two auxiliary boundary value problems involving ODEs of the form (10) but with certain adequatelychosen righthand sides , . The righthand sides and are taken to vanish at all discretization points in the original interval but not to vanish in —that is to say, the selected righthand sides correspond to sources supported in the exterior of the physical domain , see e.g. Figure 2 right. Clearly, any linear combination of the form with coincides with in the part of the spatial grid contained in . Calling and the approximate FC solutions (as described in Section 4.1) corresponding to and , it is clear that any linear combination of the form satisfies the ODE (10) with null righthand side at each one of the collocation points . It follows that, denoting by and the solutions of the linear systems
(23) 
(as shown below, the system matrix is invertible for sufficiently fine discretizations provided the functions and are chosen appropriately), the functions and are approximate FC solutions of (10) with null righthand side, which satisfy the Dirichlet boundary conditions , , and . Thus, the Fourier series
(24) 
is an approximate FC solution of ODE problem (10) that satisfies exactly the prescribed Dirichlet boundary conditions (11).
The needed invertibility of the system matrix in equation (23) for sufficiently fine discretizations indeed holds provided the functions and satisfy the conditions
(25) 
To establish this fact we let and be the exact periodic solutions of equation (10) with respective righthand sides and , and, using Lemma 1 presented in Appendix A, we show at first that the “exactsolution system matrix”, that is, the matrix that results as and are substituted by and in equation (23), is invertible. Indeed, if the exactsolution system matrix is singular then the and satisfies and for a certain . The function is then a solution of the twopoint BVP (10)–(11) with , and , and thus, it vanishes identically since, in view of Lemma 1 this problem admits a unique solution. It then follows by the smoothness and periodicity of that and . But this is a contradiction, since Lemma 1 tells us that such solutions do not exist under the hypothesis (25). Hence, the exactsolution system matrix is invertible. The invertibility for the approximatesolution system matrix (23) for sufficiently fine grids then follows from the convergence of the FC solutions and to and as the gridsize tends to zero (see Remark 6).
Remark 6.
Remark 7.
For the numerical experiments presented in this work we have used the auxiliary righthand sides
where and ; see e.g. Figure 2. As required by Lemma 1 in Appendix A, the functions and satisfy the assumption (25), and their supports, which are contained in the interval , do not intersect the interval .
Remark 8.
Note that the solution and its numerical approximation generally possess boundary layers [bender_orzag, Ch. 7] at the interval endpoints and in cases in which the coefficient function in equation (22) is small—which, in view of equations (6), (7) and (9), it certainly is for ODE problems arising from the FC PDE solver for small values of . The presence of such boundary layers can be appreciated by consideration of equation (24): the functions and provide (numerical approximations of) the boundarylayer contributions.
Remark 9.
The numerical boundarylayer FC solutions and mentioned in Remark 8 may be affected by Gibbslike ringing errors near the endpoints unless the underlying grid adequately resolves the exact boundary layers. To ensure the detection of a nonresolved boundary layer (so that a procedure can be used to guarantee that the discretization errors are not inherited by and in stiff problems), the grid size is compared to the quantities and . If the gridsize is of the order of or smaller than these quantities, the boundary layer is wellresolved by the numerical approximations and and no further action is needed. Otherwise, if is larger than and , then and do not adequately resolve the boundary layers, and an alternative treatment is necessary (Remark 16 presents details on the actual thresholds used in our numerical examples). Such an alternative method, based on use of asymptotic expansions, is presented in Section 4.3.
4.3 Boundary conditions II: Highorder asymptoticmatching expansions
As noted in the introductory paragraph of Section 4, stiff ODEs and challenging boundary layers arise in the solutions of the boundary value problem (10)–(11) as small time steps are used. In this section we present an algorithm that can resolve such boundary layers, without resort to unduly fine spatial meshes, on the basis of the method of matched asymptotic expansions [bender_orzag]. For our development of this asymptotic procedure, we introduce a small positive parameter (that is taken to equal in equation (9) and in equations (6) and (7)), we let and , and we reexpress equation (10) in the form
(26)  
(27) 
where and the variable coefficients and are bounded smooth functions defined in the interval (which satisfy the conditions for some constant ) given by
Further, using the change of variables
(28) 
equation (26) can be recast in the simpler form
(29) 
where , and
In what follows we produce, on the basis of the method of matched asymptotic expansions [bender_orzag], approximate solutions and of the homogeneous () version of equation (26), satisfying , , , ; like the functions and introduced in the previous section, and can be used to correct the boundary values of a periodic FC solution of (26). Unlike the functions introduced in the previous section, however, the asymptoticexpansion solutions and produce accurate solutions in the small regime without requiring use of fine meshes, and they can be evaluated with a fixed computational cost for arbitrarily small values of .
Remark 10.
As mentioned above in this section, the small parameter included in the variable coefficients and corresponds to the quantities and in the time discretization of the wave and heat equations, respectively. Hence, the uniform convergence of the ODE solutions as with respect to guarantees the spatial convergence of the FCAD scheme independently of the value of used in the timemarching schemes.
We lay down our procedure for evaluation of (left boundary layer); the corresponding method for (right boundary layer) is entirely analogous. Using the change of variables we define a new unknown by ; clearly satisfies equation (29) together with the boundary conditions
To obtain the approximate solution for small values of we use outer and inner solutions in the corresponding matched asymptotic expansion [bender_orzag] associated with the left boundary layer. Setting in the ODE operator and using the null boundary condition at we see that the lowest order term in the outer asymptotic expansion vanishes. In fact, as it follows from the discussion below, the outer solution vanishes to all orders:
In the inner region, in turn, we use the scaled variable and we express the inner solution and the ODE coefficient as asymptotic expansions
(30) 
In particular, the Taylor expansion of the ODE coefficient induces a corresponding expansion for the ODE operator , namely , where
It follows that the functions are solutions of the following ODE problems:
(31) 
for , and
(32) 
for . It is easy to check that , that, similarly, can be expressed in closed form in terms of adequate combinations of exponentials and polynomials (see Algorithm 1), and finally, that, as claimed above, all terms in the outer asymptotic expansion vanish.
In sum, we have
(33) 
where the functions satisfy (31) and (32); a corresponding asymptotic boundary correction function can be obtained similarly. The maximum errors resulting from these approximations satisfy and for small . Using these asymptotic boundary correction functions we then obtain an approximate FC solution of ODE problem (10), which satisfies exactly the prescribed Dirichlet boundary conditions (11), by means of the expression