Elsevier

Applied Mathematics and Computation

Volume 307, 15 August 2017, Pages 329-341
Applied Mathematics and Computation

An efficient parareal algorithm for a class of time-dependent problems with fractional Laplacian

https://doi.org/10.1016/j.amc.2017.02.012Get rights and content

Abstract

Time-dependent diffusion equations with fractional Laplacian have received considerable attention in recent years, for which numerical methods play an important role because a simple and analytic solution is often unavailable. We analyze in this paper a parareal algorithm for this kind of problem, which realizes parallel-in-time computation. The algorithm is iterative and uses the 3rd-order SDIRK (singly diagonally implicit Runge-Kutta) method with a small step-size Δt as the F-propagator and the implicit-explicit Euler method with a large step-size ΔT as the G-propagator. The two step-sizes satisfy ΔT/Δt=J with J ≥ 2 being an integer. Using the implicit-explicit Euler method as the G-propagator potentially improves the parallel efficiency, but complicates the convergence analysis. By employing some technical analysis, we provide a sharp estimate of the convergence rate, which is independent of the mesh ratio J and the distribution of the eigenvalues of the coefficient matrix. An extension of the results to problems with time-periodic conditions is also given. Several numerical experiments are carried out to verify the theoretical results.

Introduction

In this section, we introduce the time-dependent diffusion problems with fractional Laplacian studied in this paper and the parareal algorithm.

In recent years, there is a growing interest in the study of partial differential equations (PDEs) with fractional Laplacian. These equations have a wide application in many fields, such as financial mathematics [6], image processing [14], [16], turbulence [2] and many others. We consider in this paper a class of representative PDEs in this field: u(x,t)t=Δu(x,t)d(Δ)αu(x,t)+f(x,t),(x,t)Ω×(0,T),where d > 0 and α ∈ (0, 1). As adopted in [1], [3], [4], [5], [19], [20], [26], [29], [30], [37], here the symmetric spatial fractional Laplacian operator (Δ)α is defined through the eigenfunction expansion on a finite-size spatial domain (see Definition 1 in [20]). We remark that it is entirely possible to define the fractional Laplacian operator through the Fourier transform on the whole spatial domain [18], [31]. For both definitions, the idea of matrix transform proposed in [19] provides efficient spatial discretization [24], [36], [37]. The spatial discretization based on the matrix transform idea of the fractional Laplacian (Δ)α is obtained by first finding a matrix representation A of the regular (negative) Laplacian operator Δ (whether it is through finite difference, finite element, or finite volume) and then raising it to the same fractional power Aα1. Precisely, for (1.1), by introducing a mesh with m nodes and denoting the value of u(x, t) at the i-th node xi by ui(t), this treatment yields approximation u(t)(u(x1,t),u(x2,t),,u(xm,t)), where u(t) is the solution of the following ODE system u(t)+(A+dAα)u(t)=F(t),t(0,T),where F(t) contains the information of the source term f and the boundary values.

The coefficient matrix A in (1.2) can be derived, for example, through the centered finite difference scheme subjected with Neumann boundary condition (un=0 on ∂Ω), as {A=ΛΔx:=1Δx2[1112112111]m×m,1-Dcase,A=IΛΔx+ΛΔxI,2-Dcase,A=IIΛΔx+IΛΔxI+ΛΔxII,3-Dcase,where Δx denotes the step-size and m denotes the number of nodes in the spatial grid. We can also consider other boundary conditions and other spatial discretizations, such as finite element and finite volume, etc. For the chosen spatial discretization, we emphasize the following two points:

  • the matrix A is sparse and SPD (symmetric positive definite);

  • λσ(A)⇔λασ(Aα), where σ(A) denotes the spectrum of A.

Note that, in spite of the sparsity of A, the matrix Aα is dense. This implies that it would be very time-consuming to solve (1.2), if the number of spatial grids is large (i.e., the size of the matrix A is huge). Therefore, the ODE system (1.2) presents the fundamental challenges for numerical investigations of the time-dependent diffusion problems with fractional Laplacian.

To accelerate the numerical computation of (1.2), we can consider the so-called parareal algorithm proposed by Lions, Maday and Turinici [21]. The algorithm starts from splitting the whole time-interval [0, T] into N large subintervals [Tn,Tn+1] as [0,T]=n=0N1[Tn,Tn+1] with large step-size ΔT and further dividing each large subinterval [Tn,Tn+1] into J small intervals as [Tn,Tn+1]=j=0J1[Tn+j/J,Tn+(j+1)/J] with small step-size Δt. The two step-sizes satisfy ΔT/Δt=J with J ≥ 2 being an integer. Then, we choose two numerical methods, denoted by G and F, to propagate the initial solution unk on each large subinterval [Tn,Tn+1] with step-sizes ΔT and Δt, respectively. The last step is to use a novel correction scheme to improve the approximate solutions on the coarse time-grids {Tn}n=1N, as un+1k+1=G(Tn,unk+1,ΔT)+FJ(Tn,unk,Δt)G(Tn,unk,ΔT),where n=0,1,,N1, k ≥ 0 denotes the iteration index and for k=0 the initial guess {un0}n=1N are chosen randomly (or generated by a cheap numerical method). Here, G(Tn,unk,ΔT) denotes the numerical solution of (1.2) at t=Tn+1 produced by using the numerical method G with step-size ΔT and initial value unk at t=Tn. For example, if we use the Backward-Euler method as the G-propagator, G(Tn,unk,ΔT) denotes the solution u˜ of the following linear equations u˜+ΔT[(A+dAα)u˜]=unk+ΔTF(Tn+ΔT)([I+ΔT(A+dAα)]u˜=unk+ΔTF(Tn+1)).In (1.4), FJ(Tn,unk,Δt) denotes the numerical solution of (1.2) at t=Tn+1 by running J steps of the fine propagator F with small step-size Δt and initial value unk at t=Tn.

It is clear that, upon convergence we have un+1=FJ(Tn,un,Δt), i.e., the approximate solution of (1.2) at the coarse grids {Tn}n=1N will have achieved the accuracy of the fine propagator F with small step-size Δt. In the past few years, the parareal algorithm and its variants have been used in many fields, such as structural (fluid)-dynamic simulations [7], [11], optimal control [8], [23], [25], Hamiltonian computation [9], [10], [15], turbulent plasma simulations [27], [28], singularly perturbed problems [22] and time-periodic problems [13], [32], [33], etc. In spite of numerous successful applications as mentioned, the results about the convergence properties of this algorithm are less known in the literature. For the scalar linear model problem u(t)=λu(t)+f with λC, Gander and Vandewalle [12] systematically analyzed the convergence rate of this algorithm. Particularly, for any A-stable G-propagator they built a connection between the contraction factor of error of the parareal algorithm on sufficiently long time intervals and the stability functions of the propagators G and F (denoted by Rg and Rf), as supnunkunKk(z,J)supnun0un,withK(z,J):=|Rg(z)RfJ(zJ)|1|Rg(z)|andz=ΔTλ.

This result is a starting point for many subsequent research about the convergence properties of the parareal algorithm. In [25], by choosing for both the G- and F-propagators the Backward-Euler method, Mathew et al. analyzed the convergence rate of this algorithm for symmetric positive definite (SPD) system u(t)+Qu(t)=F(t) (with QRm×m being a SPD matrix) and proved the following result ρ(J):=maxzσ(ΔTQ)K(z,J)13,J2,σ(ΔTQ)[0,),where σ(ΔTQ) denotes the spectrum of the matrix ΔTQ. The quantity ρ(J) is called the convergence factor of the parareal algorithm applied to a ODE system. Since F=Backward-Euler, the accuracy of the converged solution is only of order O(Δt). Very recently, we generalized the result (1.7) to two 2nd-order F-propagators in [34], the 2nd-order SDIRK (singly diagonally implicit Runge-Kutta) method and the TR/BDF2 method (i.e., the ode23tb solver for ODEs in the software Matlab). In [35], we generalized (1.7) to a 3rd-order F-propagator, the 3rd-order SDIRK method. (For these higher order F-propagators, the proof of (1.7) requires that the mesh ratio J is an even integer.) Clearly, these results can be straightforwardly generalized to our problem (1.2), since the matrix A+dAα is SPD.

The parallel efficiency of the parareal algorithm is a different issue, but closely connected with the convergence analysis. Let Tunitf and Tunitg be respectively the computational cost of the fine and coarse propagators for moving forward one time step. Let ϵ be the anticipated tolerance. Then, as pointed out in [17, Section 5], the parallel speedup of the parareal algorithm, namely E(J)2, can be expressed as E(J):=clogρ(J)logϵ×NJJc+1+N(Jclogρ(J)logϵ,forNlarge),where c=Tunitf/Tunitg, ρ(J) denotes the convergence factor of the parareal algorithm and N=TΔT denotes the number of paralleled large subintervals. Clearly, to make E(J) large, it needs c large and ρ(J) small. For a fixed F-propagator, the requirement ‘c large’ implies that the G-propagator should be as cheap as possible. But on the other hand G should be also as stable as possible, since it is assumed to take large step-size ΔT. Therefore, previous studies mainly focus on using for the G-propagator the Backward-Euler method [7], [11], [13], [25], [34], [35]. For our problem (1.2), this choice results a linear problem with dense coefficient matrix I+ΔT(A+dAα) as shown in (1.5) for each coarse time grid.

The goal of this paper is explore the possibility for further reducing the computation time for moving forward one step of the G-propagator. Our idea lies in using the implicit explicit Euler (IMEX-Euler) method instead of the Backward-Euler method as the G-propagator. For (1.2), the IMEX-Euler method leads to the following linear system for each coarse time grid: u˜=unk+ΔTAu˜+dAαunk+F(Tn+ΔT),where the implicit and explicit computation is respectively associated with the sparse matrix A and the dense matrix Aα. Comparing (1.5) to (1.9), it is obvious that moving forward one step of the IMEX-Euler method is less expensive than the Backward-Euler method. Therefore, for a common F-propagator (such as the 3rd-order SDIRK method studied in this paper), using the IMEX-Euler method as the G-propagator possesses the potential advantage for enlarging the quantity c in (1.8). Interestingly, for α ∈ (0, 0.9), if we regard the large step-size ΔT as a free parameter and choose it appropriately, this choice also leads to a smaller convergence factor ρ(J), say ρ(J)15, for all σ(A)⊆[0, ∞) and all even integers J ≥ 8. For α ∈ [0.9, 1), using G=IMEX-Euler can give a convergence factor ρ(J)13 for all even integers J ≥ 8, which is similar to that of using G=Backward-Euler [35]. We remark that requiring the mesh ratio J to be an even integer and J ≥ 8 is a trivial restriction for real computations.

The remainder of this paper is organized as follows. In Section 2 we present some necessary notations and preparations. Section 3 presents the main result of this paper, a detailed convergence analysis of the parareal algorithm with (F,G)=(3rd-order SDIRK, IMEX-Euler). We generalize the analysis to problem (1.2) with time-periodic conditions (instead of the initial value condition). The theoretical predictions are verified in Section 5 by numerical results and we conclude this paper in Section 6.

Section snippets

Notations and preliminary preparations

For (1.2), since A is SPD we know that the coefficient matrix A+dAα is also SPD (it is diagonalisable and all the eigenvalues are positive numbers). Hence, the following result can be directly deduced from the analysis in [12, Sections 4-5].

Lemma 2.1

Letσ(ΔTA)={ΔTλ1,,ΔTλm} be the spectrum of the matrix ΔTA, Rf(z) andRg(z) be respectively the stability functions of theF- andG-propagators. If|Rg(z)|1 for all zσ(ΔTA), the error{enk:=ununk} of the parareal algorithm satisfies supnVenkρk(J)supnVen0

Convergence analysis

To analyze the convergence rate of the parareal algorithm with (F,G)=(3rd-order SDIRK, IMEX-Euler), we define y=z+γzα(z=H(y)),G(y)=1y1+H(y).Then, the function K(z,α,γ,J) given in (2.1) can be equivalently written as K(z,α,γ,J)=K^(y,α,γ,J):=|RSDIRK3J(yJ)G(y)|1|G(y)|.

Lemma 3.1

LetY0=1+(1/γ)1α. Then the function G(y) defined by (3.1) satisfies (1)y ∈ [0, Y0), G(y) > 0, G′(y) < 0; (2)G(Y0)=0,G(Y0)<0 and (3)y > Y0, G(y) < 0.

Proof

Since z=H(y) it holds that z increases as y increases and that RIMEX(z(y))=G(

Application to problems with time-periodic conditions

For problem (1.1), if we consider the time-periodic condition u(x,0)=u(x,T) instead of the initial condition u(x,0)=u0(x), we shall get the following system of ODEs {u(t)+(A+dAα)u(t)=F(t),t(0,T),u(0)=u(T)For time-periodic differential equations, Gander et al. [13] proposed the so-called PP-PC (periodic parareal algorithm with periodic coarse problem) and PP-IC (periodic parareal algorithm with initial value coarse problem) algorithms. Here, we are interested in PP-PC, because it converges

Numerical Experiments

In this section, we provide numerical results to validate the theoretical prediction of the convergence rate of the parareal algorithm using G=IMEX-Euler and F=3rd-order SDIRK. We consider the following problem with x ∈ (0, 1) and t ∈ (0, 50): {tu(x,t)=x2u(x,t)d(Δ)αu(x,t)+t(50t)cos((xx2)t0.7),u(0,t)=xu(1,t)=0.Applying the matrix transform method [3], [4], [5], [20], [29], [36] to (5.1) gives u+(A+dAα)u=t(50t)[cos((x1x12)t0.7)cos((x2x22)t0.7)cos((xmxm2)t0.7)],A=1Δx2[21121121

Conclusions

We have analyzed the convergence properties of a parareal algorithm using respectively for the G- and F-propagators the IMEX-Euler and the 3rd-order SDIRK methods, for the ODE system arising from semi-discretizing a class of time-dependent diffusion problems with fractional Laplacian. The accuracy of the converged solution is of third-order. We derived upper bound of the convergence factor, which depends on the problem parameters α and d and the large step-size ΔT, and is independent of Δt and

Acknowledgements

The author is very grateful to the anonymous referees for their careful reading of a preliminary version of the manuscript and their valuable suggestions and comments, which greatly improved the quality of this paper. This work is supported by NSFC (Nos. 11301362, 61573010 and 11671074), the Project of China Postdoctoral (2015M580777, 2016T90841), the NSF of Technology & Education of Sichuan Province (2014JQ0035, 15ZA0220) and the NSF of SUSE (2015LX01, 2016QZJ02).

References (37)

  • A. Bueno-Orovio et al.

    Fractional diffusion models of cardiac electrical propagation: role of structural heterogeneity in dispersion of repolarization

    J. R. Soc. Interface

    (2014)
  • N. Cusimano et al.

    On the order of the fractional laplacian in determining the spatio-temporal evolution of a space-fractional model of cardiac electrophysiology

    PloS one

    (2015)
  • R. Cont et al.

    A finite difference scheme for option pricing in jump diffusion and exponential lévy models

    SIAM J. Numer. Anal.

    (2005)
  • J. Cortial et al.

    A time-parallel implicit method for accelerating the solution of nonlinear structural dynamics problems

    Int. J. Numer. Meth. Engng.

    (2008)
  • X.H. Du et al.

    Inexact and truncated parareal-in-time krylov subspace methods for parabolic optimal control problems

    Electron Trans. Numer. Anal.

    (2013)
  • X. Dai et al.

    Stable parareal in time method for first- and second-order hyperbolic systems

    SIAM J. Sci. Comput.

    (2013)
  • X. Dai et al.

    Symmetric parareal algorithms for hamiltonian systems

    M2AN Math. Model Numer. Anal.

    (2012)
  • C. Farhat et al.

    Time-decomposed parallel time-integrators: theory and feasibility studies for fluid, structure, and fluid-structure applications

    Int . J. Numer. Meth. Engng.

    (2003)
  • Cited by (2)

    View full text