Spatially dispersionless, unconditionally stable FC-AD solvers for variable-coefficient PDEs

Spatially dispersionless, unconditionally stable FC-AD solvers for variable-coefficient PDEs

O. P. Bruno111Applied and Computational Mathematics, California Institute of Technology, Pasadena, CA 91125, USA – and A. Prieto222Departamento de Matemáticas, Univ. da Coruña, 15071 A Coru na, Spain.

We present fast, spatially dispersionless and unconditionally stable high-order 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, FC-based solver for two-point boundary value problems (BVP) for variable-coefficient Ordinary Differential Equations, and (iii) An Alternating Direction strategy, generalize significantly a class of FC-based solvers introduced recently for constant-coefficient PDEs. The present algorithms, which are applicable, with high-order accuracy, to variable-coefficient 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: High-order methods, Alternating Direction Implicit schemes, numerical dispersion, variable coefficient problems.

1 Introduction

We present fast, spatially dispersionless and unconditionally stable high-order 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 constant-coefficient 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, FC-based solver for two-point boundary value problems (BVP) for variable-coefficient Ordinary Differential Equations (ODE), and (iii) The Alternating Direction Implicit (ADI) methodology [peaceman55, douglas56]. One of the main enabling elements in our overall FC-AD algorithm (Fourier-Continuation Alternating-Directions) is the new solver (ii) for two-point boundary value problems with variable coefficients. Relying on preconditioners that result from inversion of oversampled finite-difference 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 FC-AD time-stepping in the variable-coefficient context. (A non-oversampled finite-difference 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 variable-coefficient FC-AD 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 high-order, essentially dispersionless character of the new FC-AD solvers by means of solutions to parabolic and hyperbolic problems. For example, the results presented in Section 7.3.1, which include FC-AD fixed-accuracy solutions at fixed numbers of points-per-wavelength for problems of sizes ranging from one to one-hundred wavelengths in size, demonstrate the spatial dispersionlessness of the FC-AD algorithm for hyperbolic equations. Further, as shown in Table 1, for example, the proposed FC-AD algorithm can evolve a solution characterized by one-million spatial unknowns in a computing time of approximately seconds per time step in a single-core run.

The remainder of this paper is organized as follows: after the introduction in Section 2.1 of the variable-coefficient 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 variable-coefficient differential equations. Sections 4 and 5 then describe our new FC-based solver for two-point BVP, which include 1) Iterative FC-based 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 FC-AD PDE solver, finally, are demonstrated in Section 7 through a variety of numerical results.

2 Preliminaries

2.1 Variable coefficients PDEs

This paper presents FC-based solvers for linear equations containing time-independent 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 variable-coefficient parabolic and hyperbolic problems




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 Peaceman-Rachford scheme [peaceman55] for diffusion problems and the non-centered scheme [lyon10] for the wave equation to the present variable-coefficient 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 hand-sides and , and the boundary values , the exact solution for of the diffusion problem (1) satisfies the Crank-Nicolson relation [peaceman55, bruno10]


Following [peaceman55], an ADI scheme can be obtained from the Crank-Nicolson iteration: denoting by the corresponding approximation of (for integer values ) and using the intermediate quantity , the ADI scheme is embodied in the equations


(with boundary condition on ) and


(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

  • Update according to

    and, finally,

  • Obtain by solving the boundary value problem


The ADI scheme thus requires solution of the one-dimensional 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 second-order local truncation error in step (D2); see [bruno10, Remark 3.2] for details. Numerical experiments we present in this paper reveal once again a second-order 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


where, once again, and . Following [lyon10] we split the stiffness term and thus obtain the ADI time-stepping scheme

  • Initialize and as

Then, for ,

  • Obtain by solving the boundary value problem,

  • Obtain as the solution of the boundary value problem

Remark 2.

Note that the expressions (D2)-(D4), and (W2)-(W3) of the alternating-direction ODEs for the heat and wave equations incorporate a division by the lowest-order variable-coefficient 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 linear-algebra 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 two-point BVP of the form


where the right-hand 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 right-hand sides of the problems (D2) and (D4), (W2) and (W3) are given by


where for (D2) and (D4)


while for (W2) and (W3)


To solve the two-point 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 Fourier-series 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 non-periodic boundary-value 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,


denotes the -dimensional space of trigonometric polynomials, and the amplitudes are given by the trapezoidal-rule expression

Note that, for conciseness, the “degree” is not explicitly displayed in the notation .

We point out that, as is well-known [hesthaven07, boyd01],

  1. An element of the set is determined uniquely by its values at the equispaced grid , where

  2. 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 non-periodic 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 left-most and right-most 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 left-most 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 FC-based PDE solvers presented in this paper rely on the FC methodology as a tool for solving two-point BVPs for ordinary differential equations with variable coefficients. As discussed in Section 5.3, to obtain an accurate solution of the relevant two-point BVPs for a given PDE problem in the present variable-coefficient context, our method requires, unlike previous FC-based 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


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


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


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 un-scaled 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 FC-based method for solution of variable-coefficient 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 non-periodic 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 non-zero 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 FC-based 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


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


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


(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.

Note that, since the coefficients and are variable, the matrix associated to the linear system (22) is not sparse. To avoid the expense associated with a direct solution of this linear system we utilize a preconditioned iterative solver, as described in Section 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 adequately-chosen right-hand sides , . The right-hand 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 right-hand 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 right-hand side at each one of the collocation points . It follows that, denoting by and the solutions of the linear systems


(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 right-hand side, which satisfy the Dirichlet boundary conditions , , and . Thus, the Fourier series


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


To establish this fact we let and be the exact -periodic solutions of equation (10) with respective right-hand sides and , and, using Lemma 1 presented in Appendix A, we show at first that the “exact-solution system matrix”, that is, the matrix that results as and are substituted by and in equation (23), is invertible. Indeed, if the exact-solution system matrix is singular then the and satisfies and for a certain . The function is then a solution of the two-point 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 exact-solution system matrix is invertible. The invertibility for the approximate-solution system matrix (23) for sufficiently fine grids then follows from the convergence of the FC solutions and to and as the grid-size tends to zero (see Remark 6).

Remark 6.

Stability and convergence proofs, which are beyond the scope of this paper, are left for future work; here we merely note that convergence of the FC solutions to exact solutions is amply demonstrated by the numerical results presented in Sections 5 and 6.

Remark 7.

For the numerical experiments presented in this work we have used the auxiliary right-hand 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 boundary-layer contributions.

Remark 9.

The numerical boundary-layer FC solutions and mentioned in Remark 8 may be affected by Gibbs-like ringing errors near the endpoints unless the underlying grid adequately resolves the exact boundary layers. To ensure the detection of a non-resolved 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 grid-size is of the order of or smaller than these quantities, the boundary layer is well-resolved 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: High-order asymptotic-matching 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 re-express equation (10) in the form


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


equation (26) can be re-cast in the simpler form


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 asymptotic-expansion 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 FC-AD scheme independently of the value of used in the time-marching 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


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:


for , and


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


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