An efficient parareal algorithm for a class of time-dependent problems with fractional Laplacian
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: 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 where u(t) is the solution of the following ODE system 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 ( on ∂Ω), as 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 as with large step-size ΔT and further dividing each large subinterval into J small intervals as with small step-size Δt. The two step-sizes satisfy with J ≥ 2 being an integer. Then, we choose two numerical methods, denoted by and to propagate the initial solution on each large subinterval 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 as where k ≥ 0 denotes the iteration index and for the initial guess are chosen randomly (or generated by a cheap numerical method). Here, denotes the numerical solution of (1.2) at produced by using the numerical method with step-size ΔT and initial value at . For example, if we use the Backward-Euler method as the -propagator, denotes the solution of the following linear equations In (1.4), denotes the numerical solution of (1.2) at by running J steps of the fine propagator with small step-size Δt and initial value at .
It is clear that, upon convergence we have i.e., the approximate solution of (1.2) at the coarse grids will have achieved the accuracy of the fine propagator 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 with Gander and Vandewalle [12] systematically analyzed the convergence rate of this algorithm. Particularly, for any A-stable -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 and (denoted by and ), as
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 - and -propagators the Backward-Euler method, Mathew et al. analyzed the convergence rate of this algorithm for symmetric positive definite (SPD) system (with being a SPD matrix) and proved the following result 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 =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 -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 -propagator, the 3rd-order SDIRK method. (For these higher order -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 is SPD.
The parallel efficiency of the parareal algorithm is a different issue, but closely connected with the convergence analysis. Let and 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 2, can be expressed as where ρ(J) denotes the convergence factor of the parareal algorithm and denotes the number of paralleled large subintervals. Clearly, to make large, it needs c large and ρ(J) small. For a fixed -propagator, the requirement ‘c large’ implies that the -propagator should be as cheap as possible. But on the other hand 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 -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 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 -propagator. Our idea lies in using the implicit explicit Euler (IMEX-Euler) method instead of the Backward-Euler method as the -propagator. For (1.2), the IMEX-Euler method leads to the following linear system for each coarse time grid: 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 -propagator (such as the 3rd-order SDIRK method studied in this paper), using the IMEX-Euler method as the -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 for all σ(A)⊆[0, ∞) and all even integers J ≥ 8. For α ∈ [0.9, 1), using =IMEX-Euler can give a convergence factor for all even integers J ≥ 8, which is similar to that of using =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 (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 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 be the spectrum of the matrix ΔTA, and be respectively the stability functions of the- and-propagators. If for all z ∈ σ(ΔTA), the error of the parareal algorithm satisfies
Convergence analysis
To analyze the convergence rate of the parareal algorithm with =(3rd-order SDIRK, IMEX-Euler), we define
Then, the function given in (2.1) can be equivalently written as
Lemma 3.1 Let. Then the function G(y) defined by (3.1) satisfies (1) ∀y ∈ [0, Y0), G(y) > 0, G′(y) < 0; (2) and (3) ∀y > Y0, G(y) < 0. Proof Since it holds that z increases as y increases and that
Application to problems with time-periodic conditions
For problem (1.1), if we consider the time-periodic condition instead of the initial condition we shall get the following system of ODEs 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 =IMEX-Euler and =3rd-order SDIRK. We consider the following problem with x ∈ (0, 1) and t ∈ (0, 50): Applying the matrix transform method [3], [4], [5], [20], [29], [36] to (5.1) gives
Conclusions
We have analyzed the convergence properties of a parareal algorithm using respectively for the - and -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)
- et al.
From diffusion-weighted MRI to anomalous diffusion imaging
Magnetic Resonance in Medicine
(2008) - et al.
A micro-macro parareal algorithm: application to singularly perturbed ordinary differential equations
SIAM J. Sci. Comput.
(2013) - et al.
Analysis of block parareal preconditioners for parabolic optimal controal problems
SIAM J. Sci. Comput.
(2010) - et al.
Mechanisms for the convergence of time-parallelized, parareal turbulent plasma simulations
J. Comput. Phys.
(2012) - et al.
A preconditioned numerical solver for stiff nonlinear reaction-diffusion equations with fractional laplacians that avoids dense matrices
J. Comput. Phys.
(2015) Krylov subspace methods for approximating functions of symmetric positive definite matrices with applications to applied statistics and anomalous diffusion
(2008)Novel analytical and numerical methods for solving fractional dynamical systems
(2005)- et al.
A concave-convex elliptic problem involving the fractional laplacian
Proc. Roy. Soc. Edinburgh Sect. A
(2013) Turbulence and dffusion. Springer Series in Synergetics
(2008)- et al.
An efficient implicit FEM scheme for fractional-in-space reaction-diffusion equations
SIAM J. Sci. Comput.
(2012)
Fractional diffusion models of cardiac electrical propagation: role of structural heterogeneity in dispersion of repolarization
J. R. Soc. Interface
On the order of the fractional laplacian in determining the spatio-temporal evolution of a space-fractional model of cardiac electrophysiology
PloS one
A finite difference scheme for option pricing in jump diffusion and exponential lévy models
SIAM J. Numer. Anal.
A time-parallel implicit method for accelerating the solution of nonlinear structural dynamics problems
Int. J. Numer. Meth. Engng.
Inexact and truncated parareal-in-time krylov subspace methods for parabolic optimal control problems
Electron Trans. Numer. Anal.
Stable parareal in time method for first- and second-order hyperbolic systems
SIAM J. Sci. Comput.
Symmetric parareal algorithms for hamiltonian systems
M2AN Math. Model Numer. Anal.
Time-decomposed parallel time-integrators: theory and feasibility studies for fluid, structure, and fluid-structure applications
Int . J. Numer. Meth. Engng.
Cited by (2)
A multigrid-reduction-in-time solver with a new two-level convergence for unsteady fractional Laplacian problems
2021, Computers and Mathematics with Applications