Method of lines transpose: High order L-stable schemes for parabolic equations using successive convolution
We present a new solver for nonlinear parabolic problems that is L-stable and achieves high order accuracy in space and time. The solver is built by first constructing a single-dimensional heat equation solver that uses fast convolution. This fundamental solver has arbitrary order of accuracy in space, and is based on the use of the Green’s function to invert a modified Helmholtz equation. Higher orders of accuracy in time are then constructed through a novel technique known as successive convolution (or resolvent expansions). These resolvent expansions facilitate our proofs of stability and convergence, and permit us to construct schemes that have provable stiff decay. The multi-dimensional solver is built by repeated application of dimensionally split independent fundamental solvers. Finally, we solve nonlinear parabolic problems by using the integrating factor method, where we apply the basic scheme to invert linear terms (that look like a heat equation), and make use of Hermite-Birkhoff interpolants to integrate the remaining nonlinear terms. Our solver is applied to several linear and nonlinear equations including heat, Allen-Cahn, and the Fitzhugh-Nagumo system of equations in one and two dimensions.
Key words. Method of lines transpose, transverse method of lines, Rothe’s method, parabolic PDEs, implicit methods, boundary integral methods, alternating direction implicit methods, ADI schemes, higher order L-stable, multiderivative
The prototypical parabolic differential equation is the heat equation. It forms a cornerstone of mathematics and physics, and its understanding is essential for defining more complicated mathematical models. Fourier introduced this equation as a means to describe transient heat flow. Fick quickly recognized its importance to particle and chemical concentrations. As a result, parabolic equations are now ubiquitous in describing diffusion processes, which are found in a vast array of physical problems, among which are reaction-diffusion models of chemical kinetics [18, 19, 34], phase field models describing morphology and pattern formation in multiphase fluids and solids [4, 3, 12], and even the volatility of stocks and bonds in mathematical finance .
Numerical solutions of (linear and nonlinear) diffusion equations have been the subject of active research for many decades. As early as the 1950’s and 60’s, it was recognized that due to the parabolic scaling, method of lines discretizations of the heat equation lead to numerically stiff systems of equations, especially for explicit time stepping. Larger time steps (on the order of the mesh spacing) can be taken with fully implicit solvers, but in practice, full matrix inversions may become difficult and costly, especially when memory is extremely limited, as was especially the case for early computers. Thus, alternate dimensionally implicit (ADI) splitting methods [13, 15, 37, 14, 16, 11], that make use of dimensional splitting and tridiagonal solvers, quickly gained popularity as part of an effort to reduce the amount of memory required to invert these systems.
Later on, memory constraints no longer defined the bottleneck for computing, and attention shifted toward methods that focused on reducing floating point operations (FLOPs), albeit with additional memory constraints. Most notable among these are Krylov methods [28, 22], boundary integral methods [20, 27], and quadrature methods [31, 23, 24, 44, 30]. However, with the advent of GPU processors, it appears that we are yet again seeing a paradigm shift towards methods that should emphasize small memory footprints, even at the expense of incurring a higher operation count. Thus, ADI-like methods, which can efficiently decompose larger problems and limit overhead communication, warrant further investigation, and these features are the motivating factor for this work.
In this paper we propose a novel numerical method for obtaining solutions to the linear heat equation, and nonlinear reaction-diffusion type equations. As an alternative to classical MOL formulations, we use the method of lines transpose (MOL), which is sometimes referred to as Rothe’s method , or the transverse method of lines . In this case, the PDE is first discretized in time, and the resulting semi-discrete (modified Helmholtz) problem can be solved using a variety of methods. From potential theory [22, 27], the solution can be constructed by discretizing boundary integrals. However, with dimensional splitting (that is related to the original ADI formulations), the MOL can be used to analytically solve simpler, one-dimensional boundary value problems, and the subsequent solution can be constructed through dimensional sweeps, resulting in an [32, 33] or [6, 7, 5] solver. Furthermore, we extend the method to higher orders of accuracy by using a novel idea referred to as successive convolution. This strategy has recently been developed in  for the wave equation by the present authors. In the present work, we not only extend the method of lines transpose to parabolic problems, but we recognize the resulting expansion as a so-called resolvent expansion [1, 21], which we leverage to prove stability and convergence of the successive convolution series. In addition, we incorporate nonlinear terms with an integrating factor method that relies on high order Hermite-Birkhoff interpolants as well as the (linear) resolvent expansions developed in this paper.
The rest of this paper is organized as follows. In Section LABEL:sec:op_calc we derive the basic scheme for the one-dimensional heat equation, which is L-stable and can achieve high orders of accuracy in space and time. In Section LABEL:sec:high-order-resolvent we describe how to obtain an arbitrary order discretization in a single dimension with resolvent expansions. In Section LABEL:sec:multiple, we describe how this can be extended to multiple dimensions, and in Section LABEL:sec:Linear we present results for linear heat in one and two dimensions. In Section LABEL:sec:Source, we describe how our approach can handle nonlinear source terms, and in Section LABEL:sec:Numerical we present numerous numerical results including Allen-Cahn and the Fitzhugh-Nagumo system of equations. Finally, we draw conclusions and point to future work in Section LABEL:sec:conclusions.
2 First order scheme in one spatial dimension
We begin by forming a semi-discrete solution to the 1D heat equation using the method of lines transpose (MOL). Let satisfy
with constant diffusion coefficient , and appropriate initial and boundary conditions. The MOL amounts to employing a finite difference scheme for the time derivative, and collocating the Laplacian term at time levels and . Following a similar approach from , we introduce a free parameter , so that the collocation has the form111In , there are a total of two time derivatives (and two space), so the right hand side depends on , and .
Next, we introduce the differential operator corresponding to the modified Helmholtz equation, defined by
After some algebra, we find that the scheme can be written as
We note that there are at least two reasonable strategies for choosing :
Maximize the order of accuracy. For example, if we choose , then the discretization is the trapezoidal rule, which is second order accurate and A-stable.
Enforce stiff decay. For example, if we choose , then the discretization is the backward Euler scheme, which is first order accurate, L-stable, yet does not maximize the order of accuracy.
Here and below, we opt for the second strategy, as the stiff decay of numerical solutions of the heat equation is of paramount importance. In Section LABEL:subsec:Stability, we develop this discussion in the context of higher order schemes that relies on a careful selection of as well as repeated applications of a single inverse operator.
Upon solving equation (LABEL:eqn:First) for , we find that the equation for the update is
that requires inverting a modified Helmholtz operator. We accomplish this analytically by using Green’s functions, from which
The coefficients and are determined by applying prescribed boundary conditions at which we describe in Section LABEL:subsec:homogenous-solution.
Alternatively, had we followed the method of lines (MOL) and first discretized (LABEL:eqn:heat) in space, then the differential operator would be replaced by an algebraic operator , and would be inverted numerically.
Although the update (LABEL:eqn:First_update) (with ) is only first order accurate, we describe in Section LABEL:sec:high-order-resolvent how to extend our procedure to arbitrary order in time.
This MOL approach has several advantages. First, the solution is now explicit, but remains unconditionally stable. Secondly, in recent work [6, 7, 5] we show that the convolution integral in Eqn. (LABEL:eqn:L_inverse) can be discretized using a fast algorithm, where is the number of spatial points. We introduce more details in Section LABEL:subsec:High_Space, wherein we update the current algorithm to achieve a user-defined accuracy of with mesh spacing . Finally, since the solution is still continuous in space, we can decouple the spatial and temporal errors, and by combining resolvent expansions with dimensional splitting, we extend the method to multiple dimensions without recoupling the errors.
Since dimensional splitting is used, all spatial quantities are computed according to a one-dimensional convolution integral of the form (LABEL:eqn:L_inverse), which is performed on a line-by-line basis, following so-called "dimensional sweeps". Since the discrete convolution is computed in complexity, the full solver scales linearly in the number of spatial points (assuming each sweep is performed in parallel).
A fully discrete scheme is obtained after a spatial discretization of (LABEL:eqn:L_inverse). The domain is partitioned into subdomains , with . The convolution operator is comprised of a particular solution, which is defined by the convolution integral
and a homogeneous solution
both of which can be constructed in operations using fast convolution. We now describe each of these in turn, starting with the first.
2.1 Spatial discretization of the particular solution
The particular solution is first split into , where
Each of these parts independently satisfy the first order "initial value problems"
where the prime denotes spatial differentiation. By symmetry, the scheme for follows from that of , which we describe. From the integrating factor method, the integral satisfies the following identity, known as exponential recursion
This expression is still exact, and indicates that only the "local" integral needs to be approximated. We therefore consider a projection of onto , the space of polynomials of degree . A local approximation
is accurate to , and defines a quadrature of the form
If standard Lagrange interpolation is used, then the polynomials can be factorized as
where , and , and is the Vandermonde matrix corresponding to the points , for . The values of and are such that , and are centered about except near the boundaries, where a one-sided stencil is required.
Substituting this factorization into (LABEL:eqn:JL) and integrating against an exponential, we find that
where the weights satisfy
If the weights are pre-computed, then the fast convolution algorithm scales as per time step, and achieves a user-defined in space. In every example shown in this work, we choose or .
2.2 Homogeneous solution
The homogeneous solution in (LABEL:eqn:homogeneous_def) is used to enforce boundary conditions. We first observe that all dependence on in the convolution integral, , in (LABEL:eqn:I_def) is on the Green’s function, which is a simple exponential function. Through direct differentiation, we obtain
Various boundary conditions at and can be enforced by solving a simple system for the unknowns and . For example, for periodic boundary conditions we assume that (at each discrete time step, )
We next enforce this assumption to hold on the scheme (LABEL:eqn:L_inverse),
where and the identities (LABEL:eqn:derivative_integral) are used to find . Solving this linear system yields
Different boundary conditions (e.g. Neumann) follow an analogous procedure that requires solving a simple linear system for and .
3 Higher order schemes from resolvent expansions
In our recent work , we apply a successive convolution approach to derive high order A-stable solvers for the wave equation. The key idea is to recognize the fact that, in view of the modified Helmholtz operator (LABEL:eqn:MHL), the second derivative can be factored as
Substitution of the second expression into (LABEL:eqn:Lap_sc) determines the second derivative completely in terms of this new operator
This shows that second order partial derivatives of a sufficiently smooth function can be approximated by truncating a resolvent expansion based on successively applying to , which is a linear combination of successive convolutions (LABEL:eqn:L_inverse).
Now, we consider a solution to the heat equation (LABEL:eqn:heat), that for simplicity we take to be infinitely smooth. We perform a Taylor expansion on , and then use the Cauchy-Kovalevskaya procedure [9, 40] to exchange temporal and spatial derivatives to yield
The term is a spatial pseudo-differential operator, and it compactly expresses the full Taylor series. Our goal is to make use of the formula (LABEL:eqn:Lap_D) to convert the Taylor series into a resolvent expansion. This can be performed term-by-term, and requires rearranging a doubly infinite sum. However, if we instead work directly with the operator defining the Taylor series, then
At first glance this expression looks quite unwieldy. However, fortune is on our side, since the generating function of the generalized Laguerre polynomials is
which bears a striking resemblance to our expansion. Indeed, if we take , substitute and , then
This expansion has been considered in the context of semigroups [1, 21], where is replaced with a general closed operator on a Hilbert space . In our notation, we restate part of Theorem 4.3 in , which is proven therein.
Let the semigroup
be approximated by
Then, for , there exists for each an integer such that for all integers , with ,
where is a constant that depends only on and .
The salient point of the theorem is that, in consideration of (c.f. Eqn. (LABEL:eqn:MHL)), the approximation error is of the form , which matches the form given by a typical Taylor method.
Finally, we truncate the resolvent expansion (LABEL:eqn:Laguerre) at order . For the heat equation, this defines the numerical method as
which has a truncation error of the form
For , these schemes (evaluated at ) are
We note that for implementation, each operator is applied successively, and is defined by
There remains one critical issue: the choice of the free parameter . In 1974, Nørsett studied a similar single-step multiderivative method for the heat equation  and he too, had a free parameter in his solver. We follow his lead on the Von-Neumann analysis based on his MOL discretization, but in this work we optimize to obtain stiff decay, whereas Nørsett chose to maximize the order of accuracy of the solver.
Consider the linear test problem
whose exact solution satisfies
Application of (LABEL:eqn:Scheme_D) to this test problem results in
The generalized Laguerre polynomials satisfy many identities, the following of which is the most pertinent:
Here, is the standard Laguerre polynomial . Following standard definitions, we say that a numerical scheme is -stable, provided in the left-half of the complex plane . Likewise, a scheme exhibits stiff decay if as . If an -stable method also exhibits stiff decay, it is -stable.
Now, observing that as , we find that
where we have used the first two expressions in (LABEL:eqn:Laguerre_Identity) to introduce a telescoping sum. We are now prepared to prove the following.
Let be an approximate solution to the heat equation (LABEL:eqn:heat), given by the successive convolution scheme (LABEL:eqn:Scheme_D). Then,
If is chosen as the smallest root of , then the scheme achieves order , but does not exhibit stiff decay.
If is chosen as the smallest root of , then the scheme achieves order , and exhibits stiff decay.
Following the first strategy, the schemes are -stable for , whereas the second strategy ensures -stability. For both strategies, -stability is achieved for , with large values of .
Proof. The proof follows by applying the maximum modulus principle coupled with (LABEL:eqn:limit). For part 1, upon examining the truncation error (LABEL:eqn:Trunc), we see that an additional order of accuracy can be gained if we choose
However, for this choice, and so stiff decay does not hold. For part 2, we instead enforce stiff decay, but then the truncation error is of order . Finally, part 3 is demonstrated by the maximum amplification factors along the imaginary axis, as shown for both strategies in Figure LABEL:fig:Stab. In particular, we observe that for .
In , the scheme was chosen to maximize the order of accuracy, implicitly leading to eliminating the first term in the truncation error (LABEL:eqn:Trunc), which is equivalent to the first strategy. However, in this work we follow the second strategy, and choose as the smallest root of to ensure stiff decay.
For comparison we record the values of chosen for each order , to those of Nørsett in Table LABEL:tab:LagZeros. For all of our solvers, we choose to be the largest possible value that still yields provable stiff decay.
|Stiff Decay||Maximal Order|
4 Resolvent expansions for multiple spatial dimensions
We extend the 1D solver to multiple spatial dimensions through the use of dimensional splitting. Our key observation is that we can use the factorization property of the exponential to perform the series expansion. For instance, in three dimensions, we have
Now, we first replace each term with the identity (LABEL:eqn:Laguerre) dimension by dimension, and then truncate the expansions which will be in terms of the univariate operators and for as defined by (LABEL:eqn:MHL), and (LABEL:eqn:Scheme_D) acting on a function . This infinite sum with three indices must then be truncated to order , and after a change of indices we find
in 3D, with the corresponding 2D operator given by
Here we adopt the notation that sums are taken over all non-negative indices that sum to , and the multinomial coefficients are defined such that and .
The proof of stability for the multi-dimensional algorithm follows directly from that of the one-dimensional case, with the same approach applied to each spatial dimension (i.e. for the 2D case, and similarly for the 3D case).
5 The heat equation
5.1 Heat equation in 1D
We first illustrate the accuracy of our method for the 1D heat equation defined in (LABEL:eqn:heat). We consider initial conditions , for with periodic boundary conditions. We integrate up to a final time of , and set . We use the fast convolution algorithm that is fourth order accurate in space (), and set the spatial grid size to be . This ensures that the dominant error in the solution is temporal. We compute errors by the -norm, and compare against the exact solution at . The result of a temporal refinement study for and is presented in Table LABEL:tab:refinement_1D_1.
5.2 Heat equation in 2D
As a second example, we present results for the 2D heat equation. We consider initial conditions , for with periodic boundary conditions. We use a uniform mesh of size . Likewise, the -error is computed by comparing against the exact solution at the final time . In Table LABEL:tab:refinement_2D_1, we present results for a temporal refinement study for orders and .
6 Reaction-diffusion systems
We next extend our method to nonlinear reaction-diffusion systems of the form
where , with , is a diffusion coefficient matrix, and the reaction term is a function of , . In the above, is a bounded domain, and we assume appropriate initial values and boundary conditions. We shall view the diffusion as being the linear part of the differential operator, and invert this linear part analytically, using successive convolution. To derive the scheme, we use operator calculus to first write
where is a pseudo-differential operator. Upon integrating (LABEL:eqn:operator_calculus) over the interval , we arrive at the update equation
where we have made use of the abbreviated notation, . On the left hand side, the diffusion terms have been collected by this pseudo-differential operator, and will be approximated using the successive convolution techniques developed above. The reaction terms on the right hand side (LABEL:eqn:first) are fully nonlinear, and we must consider nonlinear stability when choosing a method of discretization.
We first consider approximating the integral on the right hand side (LABEL:eqn:first) with the trapezoidal rule. This defines a single-step update equation, which will be second order accurate
We may also obtain a single-step third order scheme, using multiderivative integration . By replacing the integrand (LABEL:eqn:first) with a third order Hermite-Birkhoff interpolant and performing exact integration of the resulting function, we arrive at
where . The Hermite-Birkhoff interpolant that matches the integrand in (LABEL:eqn:first) at times , and , as well as its derivative at time produces the quadrature rule in (LABEL:eqn:3rd_scheme).
The proposed schemes in (LABEL:eqn:2nd_scheme) and (LABEL:eqn:3rd_scheme) produce nonlinear equations for that need to be solved at each time step. Therefore, any implicit solver will necessarily be problem dependent.
For the problems examined in this work, we make use of simple fixed-point iterative schemes. We stabilize our iterative solvers by linearizing about a background state , which depends on the problem under consideration.
6.1 A discretization of the Laplacian operator
Upon perusing the third order update equation (LABEL:eqn:3rd_scheme), we will need to use successive convolution to replace the psuedo-differential operator , as well as the Laplacian operator . This latter point has been detailed in , and so we comment briefly on it here. Using the one-dimensional expansion (LABEL:eqn:Lap_D), we observe that the two-dimensional Laplacian is similarly given by
and can be truncated at the appropriate accuracy . Here, the subscripts indicate that the convolution is only in one spatial direction, and the other variable is held fixed. Thus, is applied along horizontal lines for fixed -values, and likewise for .
7 Nonlinear numerical results
We examine in greater detail the application of our schem to the Allen-Cahn (AC) equation ,
where the reaction term is , and is a bounded domain, and satisfies homogeneous Neumann boundary conditions.
For our fixed point iteration, we linearize about the stable fixed points , which satisfy . For example, the second order scheme from (LABEL:eqn:2nd_scheme) becomes
where indicates the time step as before, and now is the iteration index. By lagging the nonlinear term , the fixed point update is made explicit. Likewise, the third order scheme from (LABEL:eqn:3rd_scheme) becomes
Here, is again understood by replacing it with a resolvent expansion, which is a truncated series of successive convolution operators.
7.2 Allen-Cahn: One-dimensional test
Here, is the speed of the traveling wave, and we choose . Additionally, we choose the delay time , so that the exact solution satisfies .
Results for a final time of are shown in Figure LABEL:fig:AC_1d, with two different time steps. The solutions agree well with the exact solution.
In Table LABEL:tab:refinement_AC, we present the -error in the numerical solution at a final time , using the exact solution (LABEL:AC_initial). We observe first order accuracy from the Backward Euler method, and the expected orders of accuracy from the second (LABEL:eqn:AC_2ndorder_iteration) and third (LABEL:eqn:AC_3rdorder_iteration) order schemes. To ensure that the temporal error is dominant, we have used the fourth order accurate scheme (eq. (LABEL:eqn:JL) with ), with to perform spatial integration in the successive convolutions.
In principle, we can achieve higher orders accuracy in space and time. The latter would require using higher order Hermite-Birkhoff interpolation to discretize the reaction term in (LABEL:eqn:first).
7.3 Allen-Cahn: Two-dimensional test
We next solve the Allen-Cahn equation in two spatial dimensions. A standard benchmark problem involves the motion of a circular interface [8, 29, 41], to which an exact solution is known in the limiting case . The radially symmetric initial conditions are defined by
which has an interfacial circle () centered at , with a radius of . This interfacial circle is unstable, and will shrink over time, as determined by the mean curvature 
Here is the velocity of the moving interface, and is the radius of the interfacial circle at time (i.e., it is the radius of the curve defined by ). By integrating (LABEL:eqn:meancurvature) with respect to time, we solve for the radius as a function of time
The moving interface problem was simulated using , , and , which are based on the parameters used in . The numerical solution is displayed in Figure LABEL:fig:AC_2d, and we observe that the interfacial circle shrinks, as is expected.
In Figure LABEL:fig:AC_2d2, we plot compare the evolution of the radii obtained by our second order scheme with the exact radius (LABEL:eqn:velocity), for two different values of the diffusion parameter . The radius is measured by taking a slice of the solution along , and then solving for the spatial point where using linear interpolation between the two closest points that satisfy and . Refinement is performed with a fixed spatial mesh , and time steps of , and . Because the radius is derived as an exact solution in the limit (i.e., ) , we observe that the smaller value of is indeed more accurate.