- Email: [email protected]
A note on efficient preconditioner of implicit Runge–Kutta methods with application to fractional diffusion equations
Applied Mathematics and Computation 351 (2019) 116–123
Contents lists available at ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
Short Communication
A note on efficient preconditioner of implicit Runge–Kutta methods with application to fractional diffusion equations Hao Chen a, Xiaoli Wang b,∗, Xiaolin Li a a b
College of Mathematics Science, Chongqing Normal University, Chongqing 401331, China College of Geography and Tourism, Chongqing Normal University, Chongqing 401331, China
a r t i c l e
i n f o
a b s t r a c t
Keywords: Fractional diffusion equations Implicit Runge–Kutta methods Iterative methods Preconditioning Kronecker product splitting
This paper is concerned with the construction of efficient preconditioning method for systems arising from implicit Runge–Kutta discretizations of time-dependent PDEs. A polynomial preconditioner is proposed based on the Kronecker product splitting iteration method introduced by Chen (BIT 54:607-721, 2014). Some useful properties on the spectrum of the preconditioned matrix are established. Numerical experiments concerning discrete solution of fractional diffusion equations are presented to show the effectiveness of the proposed polynomial preconditioner. © 2019 Published by Elsevier Inc.
1. Introduction Consider the ordinary differential equations (ODEs)
dv(t ) = −L˜v(t ) + f (t ), dt v(t0 ) = v0 ,
t ∈ [t0 , T ],
(1.1)
where v(t ), f (t ) : [t0 , T ] → RN , v0 ∈ RN , and L˜ ∈ RN×N . We are interested in the iterative solution of linear systems
Qx = b,
Q = Is IN + A L
(1.2)
that arise in the numerical integration of ODEs (1.1) based on implicit Runge–Kutta (IRK) methods (see, for instance, [17]). Here A ∈ Rs×s is the coefficient matrix of the IRK method, Is is the s × s identity matrix, L = τ L˜, τ is the integration stepsize and denotes the Kronecker product. IRK time integration schemes exhibit attractive properties such as high order of accuracy and favorable stability properties (A-stability) and they are widely used for stiff problems. However, efficient solution strategies are not easy to design. A typical s-stage IRK method will lead to coupled s × s block systems (1.2), which may be large-scale and ill-conditioned. For the numerical solution of the systems (1.2) arising from integer-order partial differential equations (PDEs), some iterative methods and preconditioning techniques have been proposed, see, for example, [1,2,4,6–9,12,19,24,28,29,31,32]. The aim of this paper is the construction of efficient preconditioner for systems (1.2) arising from IRK discretization of the fractional diffusion equations (FDEs). Since IRK methods lead to a coupled s × s block system which is in a form of the sum of two Kronecker products, the idea is to use one Kronecker product to approximate the summation of two Kronecker products. ∗
Corresponding author. E-mail addresses: [email protected] (H. Chen), [email protected] (X. Wang).
https://doi.org/10.1016/j.amc.2019.01.041 0 096-30 03/© 2019 Published by Elsevier Inc.
H. Chen, X. Wang and X. Li / Applied Mathematics and Computation 351 (2019) 116–123
117
This Kronecker product preconditioner can be derived from a Kronecker product splitting (KPS) iteration method [6]. A KPSbased polynomial preconditioner is then proposed based on that KPS iteration method and we analyze the spectrum of the corresponding preconditioned matrix. The building block of the KPS-based preconditioner has the same structure as the system matrix derived from implicit Euler discretization of the FDEs. Therefore, we can recover the (already available) high performances of the implicit Euler discretization solvers for the IRK methods. The remaining part of this paper is organized as follows. In Section 2 we review the KPS iteration method briefly. In Section 3 we present a KPS-based polynomial preconditioner based on the KPS iteration method and analyze the spectrum of the preconditioned matrix. In Section 4, numerical experiments concerning the numerical solutions of FDEs are presented to demonstrate the performance of the proposed preconditioner. Finally, some conclusions are given in Section 5. 2. Kronecker product splitting iteration In [6], an efficient iterative method, called the Kroneker product splitting iteration method, for solving linear system Qx = b was proposed. The method is based on two Kronecker product splittings of Q:
Q = (Is + α A ) IN + A (L − α IN ), and
Q = A (L + α IN ) + (Is − α A ) IN , where α is a given positive constant. Alternating between these two splittings yields the following KPS iteration
((α A + Is ) IN )xk+ 2 = (A (α IN − L ))xk + b, 1 (A (α IN + L ))xk+1 = ((α A − Is ) IN )xk+ 2 + b, 1
(2.1) 1
with k = 0, 1, . . . , and x0 arbitrary. We can eliminate the intermediate vector xk+ 2 and write the KPS iteration (2.1) in the standard form
xk+1 = G(α )xk + H (α )b, where
(2.2)
G(α ) = (α Is − A−1 )(α Is + A−1 )−1 (α IN + L )−1 (α IN − L ) ,
(2.3)
and
H (α ) = 2α (α A + Is )−1 (α IN + L )−1 . Note that G(α ) is the iteration matrix of the KPS iteration method. In fact, iteration (2.2) may also result from the splitting iteration
P (α )xk+1 = R(α )xk + b, with the splitting
Q = P (α ) − R (α ) of the coefficient matrix Q, and it is easy to show that P (α )−1 = H (α ). Straightforward computation shows that
P (α ) =
1 (α A + Is ) (α IN + L ), 2α
R (α ) =
1 (α A − Is ) (α IN − L ). 2α
(2.4)
and
And it holds that
G(α ) = P (α )−1 R(α ) = IsN − P (α )−1 Q.
(2.5)
The fixed point iteration (2.2) converges for arbitrary initial guesses and right-hand sides b to the solution x = Q −1 b if and only if ρ (G(α )) < 1, where ρ (G(α )) denotes the spectral radius of G(α ). The following results in [6] describe the convergence of the KPS iteration and the strategy to choose the optimal parameter. x0
Theorem 2.1. Assume that all the eigenvalues of A−1 have positive real parts, all the eigenvalues of L have nonnegative real parts, α is a positive constant. Then,
λ − α < 1 , ∀α > 0, ρ (G(α )) ≤ γ (α ) := max−1 λ∈σ (A ) λ + α
118
H. Chen, X. Wang and X. Li / Applied Mathematics and Computation 351 (2019) 116–123
where σ (A−1 ) denotes the spectrum of the matrix A−1 . In addition,
α ∗ := arg min α
λ − α = |λ∗ |, λ∈σ (A−1 ) λ + α max
and
γ (α ∗ ) = tan
θ∗ 2
,
where λ∗ is the eigenvalue of the matrix A−1 with the maximum argument and θ ∗ is the corresponding argument. We remark that if the underlying IRK methods are A-stable, irreducible, and have invertible coefficient matrix A, such as Gauss, RadauIIA and LobattoIIIC methods ([17]), then the real parts of the eigenvalues of A (and A−1 ) are positive. In addition, the spectrum of the matrix L that arise in a discretized FDEs typically lie in the right half of the complex plane. 3. KPS-based polynomial preconditioners Note that
Q = P (α ) − R(α ) = P (α )(IsN − G(α )). Under the assumptions of Theorem 2.1 we have that the spectral radius of G(α ) is less than one, then the inverse of Q can be written as
Q −1 = (IsN − G(α ))−1 P (α )−1 =
∞
G(α ) j P (α )−1 .
j=0
So we can take
Pk (α ) = P (α ) IsN + G(α ) + G(α )2 + · · · + G(α )k−1
−1
(3.1)
as an approximation to matrix Q. And Pk (α ) can be seen as a preconditioner to the linear system Qx = b. We call Pk (α ) the KPS-based polynomial preconditioner of matrix Q (or the KPS(k) preconditioner). Note that P1 (α ) is just the KPS preconditioner introduced in [6]. We remark that the polynomial preconditioning strategy based on splitting iteration has been studied for some different linear systems, see, for instance, [5,15,18]. At each iteration step of the KPS(k) preconditioned Krylov subspace method, one need to solve the generalized residual equation Pk (α )z = r. It follows from (3.1) that
z = IsN + G(α ) + G(α )2 + · · · + G(α )k−1 P (α )−1 r. The vector z can be obtained by the following k-step iteration:
P (α )z j = R(α )z j−1 + r, It turns out that
j = 1, 2, . . . , k.
(3.2)
zk = G(α )k z0 + IsN + G(α ) + G(α )2 + · · · + G(α )k−1 P (α )−1 r.
(3.3)
By taking z0 = 0 in (3.3) we have that zk = z. From (3.1) and (2.5) we can show that the preconditioned matrix is given by
Pk (α )−1 Q = IsN + G(α ) + G(α )2 + · · · + G(α )k−1 P (α )−1 Q = IsN − G(α )k .
(3.4)
Then by Theorem 2.1, the following result on the spectral distribution of the preconditioned matrix Pk (α )−1 Q can be easily obtained. Theorem 3.1. Let the assumptions in Theorem 2.1 be satisfied. Then, ρ (G(α )) < 1 and the eigenvalues of the preconditioned matrix Pk (α )−1 Q are located in a circle centered at (1,0) with radius (ρ (G(α )))k . We remark that a larger positive integer k in the KPS(k) preconditioner can result in a smaller radius of the spectral distribution circle around 1. However, a larger value of k leads to a higher computational cost for every iteration in a Krylov subspace method. Finding the value of k that optimizes the computational cost appears to be a difficult problem. At each iteration step of a preconditioned Krylov subspace method, one need to solve the generalized residual equation Pk (α )z = r. For computing Pk (α )−1 r, an algorithm is given in Algorithm 1. Here we use some properties of the Kronecker products, for instance,
((α A + Is ) (α IN + L ))−1 = (α A + Is )−1 (α IN + L )−1 , and
(α A + Is )−1 (α IN + L )−1 r = reshape((α IN + L )−1 R(α A + Is )−T , sN, 1 ),
where R = reshape(r, N, s ) and the operation U = reshape(V, m, n ) returns the m-by-n matrix U whose elements are taken column-wise from V.
H. Chen, X. Wang and X. Li / Applied Mathematics and Computation 351 (2019) 116–123
119
Algorithm 1 Compute z = Pk (α )−1 r. 1: 2: 3: 4: 5: 6: 7: 8: 9:
z=0 R = reshape(r, N, s ) Z = reshape(z, N, s ) For j = 1 : k Solve (α IN + L )Y = (α IN − L )Z (α A − Is )T + 2α R Solve (α A + Is )Z = Y T by the LU factorization of α A + Is Z = ZT End z = reshape(Z, sN, 1 )
4. Numerical results In this section we present numerical tests to illustrate the effectiveness of the proposed preconditioning method applied to IRK discretizations of FDEs. All experiments were performed in Matlab. We choose GMRES as our Krylov subspace method and use the Matlab function gmres as the GMRES solver. In actual computations, all runs are started from an initial vector which is chosen to be the numerical solution of the FDE in the last time step, terminated when the relative residual has been reduced by at least eight orders of magnitude (i.e., when r k 2 ≤ 10−8 r 0 2 , where rk is the residual vector after k iterations and r0 is the initial residual). Here we consider the numerical solution of the following FDEs
⎧ β β ⎪ ⎨ ∂ u(x, t ) = d+ ∂ u(x, t ) + d− ∂ u(x, t ) + f (x, t ), (x, t ) ∈ (0, 2 ) × (0, 1], ∂t ∂+ xβ ∂− xβ u(0, t ) = u(2, t ) = 0, t ∈ [0, 1], ⎪ ⎩ u(x, 0 ) = 4x2 (2 − x )2 , x ∈ [0, 2].
(4.1)
where d+ = (3 − β ), d− = 2 (3 − β ), and the source term is given by
f (x, t ) = 16π cos(4π t + −
π 2
)x2 (2 − x )2 − 32 sin(4π t +
π 2
) { x 2 −β + 2 ( 2 − x ) 2 −β
3 3 [ x 3 −β + 2 ( 2 − x ) 3 −β ] + [ x 4 −β + 2 ( 2 − x ) 4 −β ] } . 3−β (4 − β )(3 − β )
Here the left-sided ( + ) and right-sided ( − ) fractional derivatives in (4.1) are defined in Riemann–Liouville form as follows
∂ β u(x, t ) 1 ∂2 = β (2 − β ) ∂ x2 ∂+ x
∂ β u(x, t ) 1 ∂2 = (2 − β ) ∂ x2 ∂− xβ
x
u (ξ , t ) dξ , (x − ξ )β −1
2
u (ξ , t ) dξ , (ξ − x )β −1
0
x
where 1 < β < 2 is the fractional derivative order and ( · ) is the gamma function. The exact solution of this problem is u(x, t ) = 4 sin(4π t + π2 )x2 (2 − x )2 . For details of FDEs and their properties and various applications, see, e.g., [27]. By the method of lines approach, we first discretize the fractional derivative in space by the second order WSGD scheme with shifting parameters ( p, q ) = (1, 0 ) and mesh size h = 2/(N + 1 ) [30], obtaining
L˜ = h−β d+ Tβ + h−β d− TβT , where
⎡
ω1(β ) ⎢ (β ) ⎢ ω2 ⎢ . ⎢ . ⎢ . Tβ = −⎢ ⎢ .. ⎢ . ⎢ ⎢ (β ) ⎣ωN−1 ωN(β ) ω0(β ) =
β 2
(β )
g0 ,
ω0(β ) ω1(β )
ω0(β )
..
.
..
.
.
..
.
.
..
.. ..
(β ) ωN−1
ωk(β ) =
···
0 ··· .. . .. .
.
0 .. . .. . .. .
···
···
0
β 2
(β )
gk
+
ω1(β ) ω2(β )
2 − β (β ) gk−1 , 2
0
⎤
⎥ ⎥ ⎥ ⎥ ⎥ ⎥, ⎥ 0 ⎥ ⎥ ⎥ ω0(β ) ⎦ ω1(β ) 0 .. .
k ≥ 1,
120
H. Chen, X. Wang and X. Li / Applied Mathematics and Computation 351 (2019) 116–123 Table 1 Numerical results for the fifth-order RadauIIA method applied to the FDE (4.1), s = 3, τ = 0.1.
β
N+1
I Iter
CPU
Iter
CPU
Iter
CPU
Iter
CPU
Iter
CPU
Iter
CPU
1.2
27
107.2
2.62
11.8
0.56
6.0
0.49
6.0
0.66
6.2
0.68
7.8
0.73
28
256.2
21.47
11.8
1.92
6.0
1.70
6.0
2.31
6.2
2.99
7.8
3.12
29
614.1
298.61
11.8
12.61
6.0
11.04
6.0
14.50
6.2
16.56
7.8
19.76
1.5
1.8
P3 (α ∗ )
Pˆ
BTS
1467.5
2957.20
11.8
54.16
6.0
47.14
6.0
61.53
6.2
69.86
7.8
84.42
27
247.0
8.05
9.8
0.51
5.0
0.48
5.0
0.64
6.2
0.70
8.8
0.82
28
514.6
52.31
9.8
1.74
5.0
1.62
5.0
2.14
6.2
2.86
8.8
3.46
29
1235.0
631.11
9.8
11.39
5.0
10.60
5.0
13.58
6.2
17.10
8.8
22.13 94.07
210
> 20 0 0
-
9.8
48.15
5.0
45.29
5.0
59.12
6.2
70.93
8.8
27
375.0
14.54
10.0
0.52
5.0
0.46
4.8
0.63
6.2
0.71
8.6
0.81
28
759.0
88.34
10.0
1.73
5.0
1.59
4.8
2.15
6.2
2.94
8.6
3.55
29
> 20 0 0
-
10.0
11.39
5.0
10.51
4.8
13.80
6.2
16.72
8.6
21.56
210
> 20 0 0
-
10.0
49.16
5.0
45.39
4.8
58.27
6.2
70.49
8.6
92.15
(β )
gk
P2 (α ∗ )
210
and the coefficients gk (β )
P1 (α ∗ )
are the alternating fractional binomial coefficients defined as
(k − β ) k β = = (−1 ) . k (−β ) (k + 1 )
Then we use the stiffly accurate 2s − 1 order, s-stage RadauIIA method [17] to the spatial semi-discretized system. We remark that a number of preconditioning techniques have been developed for the discrete solution of FDEs, see, for instance, [3,10,11,13,16,20–23,25]. Application of the preconditioner Pk (α ) within GMRES requires solving linear subsystems with respect to coefficient matrix α A + Is and α IN + L. Since α A + Is is an s × s matrix and s is small, we solve the subsystems associated with it by means of its LU factorization (see Step 6 of Algorithm 1). We remark that matrix α IN + L has the same structure as the coefficient matrix derived from implicit Euler discretization of the FDEs. Hence, we can recover the high performances of implicit Euler discretization solvers for IRK methods. Note that α IN + L is a Toeplitz matrix in this example and (α IN + L )−1 can be expressed explicitly in the Gohberg–Semencul-type formula [14]. The matrix-vector multiplication associated with matrix (α IN + L )−1 can be obtained by the FFT (see [21] for details). We remark that the matrix-vector product associated with the Toeplitz matrix L can be obtained by the FFT (see [26]). For numerical comparison, we also present results obtained by applying a block triangular splitting (BTS) method proposed in [1]. The idea is to transform the linear system Qx = b to
Qˆ y := (Is IN + Aˆ L )y = bˆ ,
x = (V −1 IN )y,
where Aˆ = VAV −1 and bˆ = (V IN )b. Then we choose the transformation matrix V (see [1] for details on obtaining V) such that Aˆ = LˆUˆ , where Uˆ is upper triangular with unit diagonal entries and Lˆ is lower triangular with the elements on the main diagonal all equal to l = s det(A ). Then the BTS iteration is defined as
Pˆy(k+1) = (Pˆ − Qˆ )y(k ) + bˆ = ((Lˆ − Aˆ ) L )y(k ) + bˆ ,
k = 0, 1, . . . ,
(4.2)
where
Pˆ = Is IN + Lˆ L.
(4.3)
Note that we can rewrite the BTS iteration (4.2) in correction form:
y(k+1) = y(k ) + Pˆ−1 r (k ) ,
r (k ) = bˆ − Qˆ y(k ) ,
(4.4)
and the splitting matrix Pˆ can be seen as a preconditioner for matrix Qˆ . In this section we present results obtained by applying the BTS method both as a stationary iteration scheme and as a preconditioner for GMRES. We remark that both the computation of the BTS iteration (4.4) and the application of the BTS preconditioner Pˆ within GMRES require solving linear subsystems with respect to coefficient matrix IN + lL (the diagonal block of Pˆ). The matrix-vector multiplication associated with (IN + lL )−1 is also expressed in the Gohberg-Semencul-type formula [14] and then computed by the FFT. In Table 1 we present the average number of iterations (Iter) per time step and the total CPU time (CPU) measured in seconds for solving the FDE (4.1). We compare the performances of the GMRES without preconditioner (I), GMRES with the KPS-based polynomial preconditioners Pk (α ∗ ) for k = 1, 2, 3, GMRES with the BTS preconditioner Pˆ and the stationary BTS iteration. We can see that all the KPS-based polynomial preconditioners Pk (α ∗ ) can largely reduce the number of iteration
H. Chen, X. Wang and X. Li / Applied Mathematics and Computation 351 (2019) 116–123
121
Fig. 4.1. The eigenvalue distribution of the matrices Q, P1 (α ∗ )−1 Q, P2 (α ∗ )−1 Q and P3 (α ∗ )−1 Q for the fifth-order RadauIIA discretization of the FDE (4.1) (β = 1.8, s = 3, N = 28 − 1, τ = 0.1).
steps of the preconditioned GMRES method. The preconditioner P2 (α ∗ ) performs better than the preconditioners P1 (α ∗ ) and P3 (α ∗ ) in terms of CPU times. The numerical results show that the BTS-preconditioned GMRES performs better than the stationary BTS iteration. Furthermore, we note that the performance of the KPS-based polynomial preconditioners Pk (α ∗ ) (k = 1, 2, 3) is better than that of the BTS preconditioner tested in terms of computational times. In Fig. 4.1 we display the eigenvalues of the matrix Q before and after preconditioning with Pk (α ∗ ) (k = 1, 2, 3). These figures clearly show that the matrix without preconditioning is very ill conditioned while the matrices with the proposed preconditioning have tightly clustered eigenvalues and, thus, are well conditioned. In Fig. 4.2 we summarize results of the 2s − 1-order RadauIIA methods (RadauIIA-2s − 1) with s = 1, 2, 3, 4 for the FDE (4.1) with β = 1.5. Note that the first-order RadauIIA method is just the implicit Euler method. The spatial fractional operators of the FDE are discretized by WSGD scheme with mesh size h = 2−9 (N = 210 − 1). All the linear systems are solved by GMRES with preconditioner P2 (α ∗ ) (or Pˆ) except the ones associated with implicit Euler method, which are solved directly by the Gohberg-Semencul-type formula. We compare the accuracy against time step size and the total CPU time. Here, the error is measured in the discrete L2 -norm at t = 1. It is shown that higher-order RadauIIA methods require fewer steps to achieve a given accuracy. They are also faster methods provided that the given accuracy is high enough. From Fig. 4.2, we also see that the performance of the BTS preconditioner is better than the KPS-based preconditioner P2 (α ∗ ) in terms of computational cost for RadauIIA-3 method, while the preconditioner P2 (α ∗ ) performs better than the BTS preconditioner for RadauIIA-5 and RadauIIA-7 methods.
122
H. Chen, X. Wang and X. Li / Applied Mathematics and Computation 351 (2019) 116–123
Fig. 4.2. Accuracy versus time-step and CPU time for four RadauIIA methods applied to the FDE (4.1) with β = 1.5.
5. Conclusions In this paper we have proposed a KPS-based polynomial preconditioning strategy for systems which arise from implicit Runge-Kutta discretizations of the FDEs. The proposed preconditioner is constructed based on the KPS iteration method. One component of the KPS-based polynomial preconditioner has the same structure as the system derived from implicit Euler discretization of the FDEs. So we can reuse the high performance of implicit Euler discretization solver as the building block for our preconditioner. Theoretical and numerical results show that the proposed preconditioning method is effective and efficient. Acknowledgments The authors are very grateful to the referees for their constructive comments and valuable suggestions, which greatly improved the original manuscript of the paper. This work was supported by the National Natural Science Foundation of China (grant nos. 11301575 and 11471063), the Scientific and Technological Research Program of Chongqing Municipal Education Commission (grant nos. KJ1703055 and KJZD-M201800501), the Natural Science Foundation Project of CQ CSTC (nos. cstc2018jcyjAX0113 and cstc2018jcyjAX0794), and the Talent Project of Chongqing Normal University (grant no. 02030307– 0054). References [1] P. Amodio, L. Brugnano, A note on the efficient implementation of implicit methods for ODEs, J. Comput. Appl. Math. 87 (1997) 1–9. [2] O. Axelsson, R. Blaheta, R. Kohut, Preconditioning methods for high-order strongly stable time integration methods with an application for a DAE problem, Numer. Linear Algebra. Appl. 22 (2015) 930–949. [3] Z.Z. Bai, K.Y. Lu, J.Y. Pan, Diagonal and toeplitz splitting iteration methods for diagonal-plus-toeplitz linear systems from spatial fractional diffusion equations, Numer. Linear Algebra. Appl. (2017), doi:10.1002/nla.2093. [4] L. Brugnano, F. Iavernaro, C. Magherini, Efficient implementation of radau collocation methods, Appl. Numer. Math. 87 (2015) 100–113. [5] Z.H. Cao, Y. Wang, On validity of m-step multisplitting preconditioners for linear systems, Appl. Math. Comput. 126 (2002) 199–211. [6] H. Chen, A splitting preconditioner for the iterative solution of implicit Runge-Kutta and boundary value methods, BIT 54 (2014) 607–621. [7] H. Chen, Generalized kronecker product splitting iteration for the solution of implicit Runge-Kutta and boundary value methods, Numer. Linear Algebra. Appl. 22 (2015) 357–370. [8] H. Chen, Kronecker product splitting preconditioners for implicit runge-kutta discretizations of viscous wave equations, Appl. Math. Model. 40 (2016) 4429–4440. [9] H. Chen, A splitting preconditioner for implicit runge-kutta discretizations of a partial differential-algebraic equation, Numer. Algor. 73 (2016) 1037–1054. [10] H. Chen, W. Lv, T.T. Zhang, A kronecker product splitting preconditioner for two-dimensional space-fractional diffusion equations, J. Comput. Phys. 360 (2018) 1–14. [11] H. Chen, T.T. Zhang, W. Lv, Block preconditioning strategies for time-space fractional diffusion equations, Appl. Math. Comput. 337 (2018) 41–53. [12] G.J. Cooper, J.C. Butcher, An iteration scheme for implicit Runge-Kutta methods, IMA J. Numer. Anal. 3 (1983) 127–140. [13] M. Donatelli, M. Mazza, S. Serra-Capizzano, Spectral analysis and structure preserving preconditioners for fractional diffusion equations, J. Comput. Phys. 307 (2016) 262–279. [14] I. Gohberg, V. Olshevsky, Circulants, displacements and decompositions of matrices, Integral Equ. Oper. Theory 15 (1992) 730–743. [15] X.M. Gu, T.Z. Huang, H.B. Li, L. Li, W.H. Luo, On k-step CSCS-based polynomial preconditioners for toeplitz linear systems with application to fractional diffusion equations, Appl. Math. Lett. 42 (2015a) 53–58. [16] X.M. Gu, T.Z. Huang, X.L. Zhao, H.B. Li, L. Li, Strang-type preconditioners for solving fractional diffusion equations by boundary value methods, J. Comput. Appl. Math. 277 (2015b) 73–86. [17] E. Hairer, G. Wanner, Solving Ordinary Differential Equations II. Stiff and Differential Algebraic Problems, Springer-Verlag, Berlin, 1996. [18] Y.M. Huang, On m-step hermitian and skew-hermitian splitting preconditioning methods, J. Eng. Math. 93 (2015) 77–86. [19] L.O. Jay, T. Braconnier, A parallelizable preconditioner for the iterative solution of implicit Runge-Kutta type methods, J. Comput. Appl. Math. 111 (1999) 63–76. [20] S.L. Lei, H.W. Sun, A circulant preconditioner for fractional diffusion equations, J. Comput. Phys. 242 (2013) 715–725.
H. Chen, X. Wang and X. Li / Applied Mathematics and Computation 351 (2019) 116–123
123
[21] S.L. Lei, Y.C. Huang, Fast algorithms for high-order numerical methods for space-fractional diffusion equations, Int. J. Comput. Math. 94 (2017) 1062–1078. [22] X.L. Lin, M.K. Ng, H.W. Sun, A splitting preconditioner for toeplitz-like linear systems arising from fractional diffusion equations, SIAM J. Matrix Anal. Appl. 38 (2017) 1580–1614. [23] X.L. Lin, M.K. Ng, H.W. Sun, Efficient preconditioner of one-sided space fractional diffusion equation, BIT 58 (2018) 729–748. [24] K.A. Mardal, T.K. Nilssen, G.A. Staff, Order optimal preconditioners for implicit runge-kutta schemes applied to parabolic PDEs, SIAM J. Sci. Comput. 29 (2007) 361–375. [25] H. Moghaderi, M. Dehghan, M. Donatelli, M. Mazza, Spectral analysis and multigrid preconditioners for two-dimensional space-fractional diffusion equations, J. Comput. Phys. 350 (2017) 992–1011. [26] M. Ng, Iterative methods for Toeplitz systems, Oxford University Press, 2004. [27] I. Podlubny, Fractional differential equations, Math. Sci. Eng. 198 (1999). Academic Press, San Diego. [28] S. Perez-Rodriguez, S. Gonzalez-Pinto, B.P. Sommeijer, An iterated radau method for time-dependent PDEs, J. Comput. Appl. Math. 231 (2009) 49–66. [29] G.A. Staff, K.A. Mardal, T.K. Nilssen, Preconditioning of full implicit Runge-Kutta schemes for parabolic PDEs, Model. Identif. Control 27 (2006) 109–123. [30] W. Tian, H. Zhou, W. Deng, A class of second order difference approximation for solving space-fractional diffusion equations, Math. Comp. 84 (2015) 1703–1727. [31] P.J. van der Houwen, J.J.B. de Swart, Triangularly implicit iteration methods for ODE-IVP solvers, SIAM J. Sci. Comput. 18 (1997) 41–55. [32] J.V. Lent, S. Vandewalle, Multigrid methods for implicit runge-kutta and boundary value method discretizations of PDEs, SIAM J. Sci. Comput. 27 (2005) 67–92.
Contents lists available at ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
Short Communication
A note on efficient preconditioner of implicit Runge–Kutta methods with application to fractional diffusion equations Hao Chen a, Xiaoli Wang b,∗, Xiaolin Li a a b
College of Mathematics Science, Chongqing Normal University, Chongqing 401331, China College of Geography and Tourism, Chongqing Normal University, Chongqing 401331, China
a r t i c l e
i n f o
a b s t r a c t
Keywords: Fractional diffusion equations Implicit Runge–Kutta methods Iterative methods Preconditioning Kronecker product splitting
This paper is concerned with the construction of efficient preconditioning method for systems arising from implicit Runge–Kutta discretizations of time-dependent PDEs. A polynomial preconditioner is proposed based on the Kronecker product splitting iteration method introduced by Chen (BIT 54:607-721, 2014). Some useful properties on the spectrum of the preconditioned matrix are established. Numerical experiments concerning discrete solution of fractional diffusion equations are presented to show the effectiveness of the proposed polynomial preconditioner. © 2019 Published by Elsevier Inc.
1. Introduction Consider the ordinary differential equations (ODEs)
dv(t ) = −L˜v(t ) + f (t ), dt v(t0 ) = v0 ,
t ∈ [t0 , T ],
(1.1)
where v(t ), f (t ) : [t0 , T ] → RN , v0 ∈ RN , and L˜ ∈ RN×N . We are interested in the iterative solution of linear systems
Qx = b,
Q = Is IN + A L
(1.2)
that arise in the numerical integration of ODEs (1.1) based on implicit Runge–Kutta (IRK) methods (see, for instance, [17]). Here A ∈ Rs×s is the coefficient matrix of the IRK method, Is is the s × s identity matrix, L = τ L˜, τ is the integration stepsize and denotes the Kronecker product. IRK time integration schemes exhibit attractive properties such as high order of accuracy and favorable stability properties (A-stability) and they are widely used for stiff problems. However, efficient solution strategies are not easy to design. A typical s-stage IRK method will lead to coupled s × s block systems (1.2), which may be large-scale and ill-conditioned. For the numerical solution of the systems (1.2) arising from integer-order partial differential equations (PDEs), some iterative methods and preconditioning techniques have been proposed, see, for example, [1,2,4,6–9,12,19,24,28,29,31,32]. The aim of this paper is the construction of efficient preconditioner for systems (1.2) arising from IRK discretization of the fractional diffusion equations (FDEs). Since IRK methods lead to a coupled s × s block system which is in a form of the sum of two Kronecker products, the idea is to use one Kronecker product to approximate the summation of two Kronecker products. ∗
Corresponding author. E-mail addresses: [email protected] (H. Chen), [email protected] (X. Wang).
https://doi.org/10.1016/j.amc.2019.01.041 0 096-30 03/© 2019 Published by Elsevier Inc.
H. Chen, X. Wang and X. Li / Applied Mathematics and Computation 351 (2019) 116–123
117
This Kronecker product preconditioner can be derived from a Kronecker product splitting (KPS) iteration method [6]. A KPSbased polynomial preconditioner is then proposed based on that KPS iteration method and we analyze the spectrum of the corresponding preconditioned matrix. The building block of the KPS-based preconditioner has the same structure as the system matrix derived from implicit Euler discretization of the FDEs. Therefore, we can recover the (already available) high performances of the implicit Euler discretization solvers for the IRK methods. The remaining part of this paper is organized as follows. In Section 2 we review the KPS iteration method briefly. In Section 3 we present a KPS-based polynomial preconditioner based on the KPS iteration method and analyze the spectrum of the preconditioned matrix. In Section 4, numerical experiments concerning the numerical solutions of FDEs are presented to demonstrate the performance of the proposed preconditioner. Finally, some conclusions are given in Section 5. 2. Kronecker product splitting iteration In [6], an efficient iterative method, called the Kroneker product splitting iteration method, for solving linear system Qx = b was proposed. The method is based on two Kronecker product splittings of Q:
Q = (Is + α A ) IN + A (L − α IN ), and
Q = A (L + α IN ) + (Is − α A ) IN , where α is a given positive constant. Alternating between these two splittings yields the following KPS iteration
((α A + Is ) IN )xk+ 2 = (A (α IN − L ))xk + b, 1 (A (α IN + L ))xk+1 = ((α A − Is ) IN )xk+ 2 + b, 1
(2.1) 1
with k = 0, 1, . . . , and x0 arbitrary. We can eliminate the intermediate vector xk+ 2 and write the KPS iteration (2.1) in the standard form
xk+1 = G(α )xk + H (α )b, where
(2.2)
G(α ) = (α Is − A−1 )(α Is + A−1 )−1 (α IN + L )−1 (α IN − L ) ,
(2.3)
and
H (α ) = 2α (α A + Is )−1 (α IN + L )−1 . Note that G(α ) is the iteration matrix of the KPS iteration method. In fact, iteration (2.2) may also result from the splitting iteration
P (α )xk+1 = R(α )xk + b, with the splitting
Q = P (α ) − R (α ) of the coefficient matrix Q, and it is easy to show that P (α )−1 = H (α ). Straightforward computation shows that
P (α ) =
1 (α A + Is ) (α IN + L ), 2α
R (α ) =
1 (α A − Is ) (α IN − L ). 2α
(2.4)
and
And it holds that
G(α ) = P (α )−1 R(α ) = IsN − P (α )−1 Q.
(2.5)
The fixed point iteration (2.2) converges for arbitrary initial guesses and right-hand sides b to the solution x = Q −1 b if and only if ρ (G(α )) < 1, where ρ (G(α )) denotes the spectral radius of G(α ). The following results in [6] describe the convergence of the KPS iteration and the strategy to choose the optimal parameter. x0
Theorem 2.1. Assume that all the eigenvalues of A−1 have positive real parts, all the eigenvalues of L have nonnegative real parts, α is a positive constant. Then,
λ − α < 1 , ∀α > 0, ρ (G(α )) ≤ γ (α ) := max−1 λ∈σ (A ) λ + α
118
H. Chen, X. Wang and X. Li / Applied Mathematics and Computation 351 (2019) 116–123
where σ (A−1 ) denotes the spectrum of the matrix A−1 . In addition,
α ∗ := arg min α
λ − α = |λ∗ |, λ∈σ (A−1 ) λ + α max
and
γ (α ∗ ) = tan
θ∗ 2
,
where λ∗ is the eigenvalue of the matrix A−1 with the maximum argument and θ ∗ is the corresponding argument. We remark that if the underlying IRK methods are A-stable, irreducible, and have invertible coefficient matrix A, such as Gauss, RadauIIA and LobattoIIIC methods ([17]), then the real parts of the eigenvalues of A (and A−1 ) are positive. In addition, the spectrum of the matrix L that arise in a discretized FDEs typically lie in the right half of the complex plane. 3. KPS-based polynomial preconditioners Note that
Q = P (α ) − R(α ) = P (α )(IsN − G(α )). Under the assumptions of Theorem 2.1 we have that the spectral radius of G(α ) is less than one, then the inverse of Q can be written as
Q −1 = (IsN − G(α ))−1 P (α )−1 =
∞
G(α ) j P (α )−1 .
j=0
So we can take
Pk (α ) = P (α ) IsN + G(α ) + G(α )2 + · · · + G(α )k−1
−1
(3.1)
as an approximation to matrix Q. And Pk (α ) can be seen as a preconditioner to the linear system Qx = b. We call Pk (α ) the KPS-based polynomial preconditioner of matrix Q (or the KPS(k) preconditioner). Note that P1 (α ) is just the KPS preconditioner introduced in [6]. We remark that the polynomial preconditioning strategy based on splitting iteration has been studied for some different linear systems, see, for instance, [5,15,18]. At each iteration step of the KPS(k) preconditioned Krylov subspace method, one need to solve the generalized residual equation Pk (α )z = r. It follows from (3.1) that
z = IsN + G(α ) + G(α )2 + · · · + G(α )k−1 P (α )−1 r. The vector z can be obtained by the following k-step iteration:
P (α )z j = R(α )z j−1 + r, It turns out that
j = 1, 2, . . . , k.
(3.2)
zk = G(α )k z0 + IsN + G(α ) + G(α )2 + · · · + G(α )k−1 P (α )−1 r.
(3.3)
By taking z0 = 0 in (3.3) we have that zk = z. From (3.1) and (2.5) we can show that the preconditioned matrix is given by
Pk (α )−1 Q = IsN + G(α ) + G(α )2 + · · · + G(α )k−1 P (α )−1 Q = IsN − G(α )k .
(3.4)
Then by Theorem 2.1, the following result on the spectral distribution of the preconditioned matrix Pk (α )−1 Q can be easily obtained. Theorem 3.1. Let the assumptions in Theorem 2.1 be satisfied. Then, ρ (G(α )) < 1 and the eigenvalues of the preconditioned matrix Pk (α )−1 Q are located in a circle centered at (1,0) with radius (ρ (G(α )))k . We remark that a larger positive integer k in the KPS(k) preconditioner can result in a smaller radius of the spectral distribution circle around 1. However, a larger value of k leads to a higher computational cost for every iteration in a Krylov subspace method. Finding the value of k that optimizes the computational cost appears to be a difficult problem. At each iteration step of a preconditioned Krylov subspace method, one need to solve the generalized residual equation Pk (α )z = r. For computing Pk (α )−1 r, an algorithm is given in Algorithm 1. Here we use some properties of the Kronecker products, for instance,
((α A + Is ) (α IN + L ))−1 = (α A + Is )−1 (α IN + L )−1 , and
(α A + Is )−1 (α IN + L )−1 r = reshape((α IN + L )−1 R(α A + Is )−T , sN, 1 ),
where R = reshape(r, N, s ) and the operation U = reshape(V, m, n ) returns the m-by-n matrix U whose elements are taken column-wise from V.
H. Chen, X. Wang and X. Li / Applied Mathematics and Computation 351 (2019) 116–123
119
Algorithm 1 Compute z = Pk (α )−1 r. 1: 2: 3: 4: 5: 6: 7: 8: 9:
z=0 R = reshape(r, N, s ) Z = reshape(z, N, s ) For j = 1 : k Solve (α IN + L )Y = (α IN − L )Z (α A − Is )T + 2α R Solve (α A + Is )Z = Y T by the LU factorization of α A + Is Z = ZT End z = reshape(Z, sN, 1 )
4. Numerical results In this section we present numerical tests to illustrate the effectiveness of the proposed preconditioning method applied to IRK discretizations of FDEs. All experiments were performed in Matlab. We choose GMRES as our Krylov subspace method and use the Matlab function gmres as the GMRES solver. In actual computations, all runs are started from an initial vector which is chosen to be the numerical solution of the FDE in the last time step, terminated when the relative residual has been reduced by at least eight orders of magnitude (i.e., when r k 2 ≤ 10−8 r 0 2 , where rk is the residual vector after k iterations and r0 is the initial residual). Here we consider the numerical solution of the following FDEs
⎧ β β ⎪ ⎨ ∂ u(x, t ) = d+ ∂ u(x, t ) + d− ∂ u(x, t ) + f (x, t ), (x, t ) ∈ (0, 2 ) × (0, 1], ∂t ∂+ xβ ∂− xβ u(0, t ) = u(2, t ) = 0, t ∈ [0, 1], ⎪ ⎩ u(x, 0 ) = 4x2 (2 − x )2 , x ∈ [0, 2].
(4.1)
where d+ = (3 − β ), d− = 2 (3 − β ), and the source term is given by
f (x, t ) = 16π cos(4π t + −
π 2
)x2 (2 − x )2 − 32 sin(4π t +
π 2
) { x 2 −β + 2 ( 2 − x ) 2 −β
3 3 [ x 3 −β + 2 ( 2 − x ) 3 −β ] + [ x 4 −β + 2 ( 2 − x ) 4 −β ] } . 3−β (4 − β )(3 − β )
Here the left-sided ( + ) and right-sided ( − ) fractional derivatives in (4.1) are defined in Riemann–Liouville form as follows
∂ β u(x, t ) 1 ∂2 = β (2 − β ) ∂ x2 ∂+ x
∂ β u(x, t ) 1 ∂2 = (2 − β ) ∂ x2 ∂− xβ
x
u (ξ , t ) dξ , (x − ξ )β −1
2
u (ξ , t ) dξ , (ξ − x )β −1
0
x
where 1 < β < 2 is the fractional derivative order and ( · ) is the gamma function. The exact solution of this problem is u(x, t ) = 4 sin(4π t + π2 )x2 (2 − x )2 . For details of FDEs and their properties and various applications, see, e.g., [27]. By the method of lines approach, we first discretize the fractional derivative in space by the second order WSGD scheme with shifting parameters ( p, q ) = (1, 0 ) and mesh size h = 2/(N + 1 ) [30], obtaining
L˜ = h−β d+ Tβ + h−β d− TβT , where
⎡
ω1(β ) ⎢ (β ) ⎢ ω2 ⎢ . ⎢ . ⎢ . Tβ = −⎢ ⎢ .. ⎢ . ⎢ ⎢ (β ) ⎣ωN−1 ωN(β ) ω0(β ) =
β 2
(β )
g0 ,
ω0(β ) ω1(β )
ω0(β )
..
.
..
.
.
..
.
.
..
.. ..
(β ) ωN−1
ωk(β ) =
···
0 ··· .. . .. .
.
0 .. . .. . .. .
···
···
0
β 2
(β )
gk
+
ω1(β ) ω2(β )
2 − β (β ) gk−1 , 2
0
⎤
⎥ ⎥ ⎥ ⎥ ⎥ ⎥, ⎥ 0 ⎥ ⎥ ⎥ ω0(β ) ⎦ ω1(β ) 0 .. .
k ≥ 1,
120
H. Chen, X. Wang and X. Li / Applied Mathematics and Computation 351 (2019) 116–123 Table 1 Numerical results for the fifth-order RadauIIA method applied to the FDE (4.1), s = 3, τ = 0.1.
β
N+1
I Iter
CPU
Iter
CPU
Iter
CPU
Iter
CPU
Iter
CPU
Iter
CPU
1.2
27
107.2
2.62
11.8
0.56
6.0
0.49
6.0
0.66
6.2
0.68
7.8
0.73
28
256.2
21.47
11.8
1.92
6.0
1.70
6.0
2.31
6.2
2.99
7.8
3.12
29
614.1
298.61
11.8
12.61
6.0
11.04
6.0
14.50
6.2
16.56
7.8
19.76
1.5
1.8
P3 (α ∗ )
Pˆ
BTS
1467.5
2957.20
11.8
54.16
6.0
47.14
6.0
61.53
6.2
69.86
7.8
84.42
27
247.0
8.05
9.8
0.51
5.0
0.48
5.0
0.64
6.2
0.70
8.8
0.82
28
514.6
52.31
9.8
1.74
5.0
1.62
5.0
2.14
6.2
2.86
8.8
3.46
29
1235.0
631.11
9.8
11.39
5.0
10.60
5.0
13.58
6.2
17.10
8.8
22.13 94.07
210
> 20 0 0
-
9.8
48.15
5.0
45.29
5.0
59.12
6.2
70.93
8.8
27
375.0
14.54
10.0
0.52
5.0
0.46
4.8
0.63
6.2
0.71
8.6
0.81
28
759.0
88.34
10.0
1.73
5.0
1.59
4.8
2.15
6.2
2.94
8.6
3.55
29
> 20 0 0
-
10.0
11.39
5.0
10.51
4.8
13.80
6.2
16.72
8.6
21.56
210
> 20 0 0
-
10.0
49.16
5.0
45.39
4.8
58.27
6.2
70.49
8.6
92.15
(β )
gk
P2 (α ∗ )
210
and the coefficients gk (β )
P1 (α ∗ )
are the alternating fractional binomial coefficients defined as
(k − β ) k β = = (−1 ) . k (−β ) (k + 1 )
Then we use the stiffly accurate 2s − 1 order, s-stage RadauIIA method [17] to the spatial semi-discretized system. We remark that a number of preconditioning techniques have been developed for the discrete solution of FDEs, see, for instance, [3,10,11,13,16,20–23,25]. Application of the preconditioner Pk (α ) within GMRES requires solving linear subsystems with respect to coefficient matrix α A + Is and α IN + L. Since α A + Is is an s × s matrix and s is small, we solve the subsystems associated with it by means of its LU factorization (see Step 6 of Algorithm 1). We remark that matrix α IN + L has the same structure as the coefficient matrix derived from implicit Euler discretization of the FDEs. Hence, we can recover the high performances of implicit Euler discretization solvers for IRK methods. Note that α IN + L is a Toeplitz matrix in this example and (α IN + L )−1 can be expressed explicitly in the Gohberg–Semencul-type formula [14]. The matrix-vector multiplication associated with matrix (α IN + L )−1 can be obtained by the FFT (see [21] for details). We remark that the matrix-vector product associated with the Toeplitz matrix L can be obtained by the FFT (see [26]). For numerical comparison, we also present results obtained by applying a block triangular splitting (BTS) method proposed in [1]. The idea is to transform the linear system Qx = b to
Qˆ y := (Is IN + Aˆ L )y = bˆ ,
x = (V −1 IN )y,
where Aˆ = VAV −1 and bˆ = (V IN )b. Then we choose the transformation matrix V (see [1] for details on obtaining V) such that Aˆ = LˆUˆ , where Uˆ is upper triangular with unit diagonal entries and Lˆ is lower triangular with the elements on the main diagonal all equal to l = s det(A ). Then the BTS iteration is defined as
Pˆy(k+1) = (Pˆ − Qˆ )y(k ) + bˆ = ((Lˆ − Aˆ ) L )y(k ) + bˆ ,
k = 0, 1, . . . ,
(4.2)
where
Pˆ = Is IN + Lˆ L.
(4.3)
Note that we can rewrite the BTS iteration (4.2) in correction form:
y(k+1) = y(k ) + Pˆ−1 r (k ) ,
r (k ) = bˆ − Qˆ y(k ) ,
(4.4)
and the splitting matrix Pˆ can be seen as a preconditioner for matrix Qˆ . In this section we present results obtained by applying the BTS method both as a stationary iteration scheme and as a preconditioner for GMRES. We remark that both the computation of the BTS iteration (4.4) and the application of the BTS preconditioner Pˆ within GMRES require solving linear subsystems with respect to coefficient matrix IN + lL (the diagonal block of Pˆ). The matrix-vector multiplication associated with (IN + lL )−1 is also expressed in the Gohberg-Semencul-type formula [14] and then computed by the FFT. In Table 1 we present the average number of iterations (Iter) per time step and the total CPU time (CPU) measured in seconds for solving the FDE (4.1). We compare the performances of the GMRES without preconditioner (I), GMRES with the KPS-based polynomial preconditioners Pk (α ∗ ) for k = 1, 2, 3, GMRES with the BTS preconditioner Pˆ and the stationary BTS iteration. We can see that all the KPS-based polynomial preconditioners Pk (α ∗ ) can largely reduce the number of iteration
H. Chen, X. Wang and X. Li / Applied Mathematics and Computation 351 (2019) 116–123
121
Fig. 4.1. The eigenvalue distribution of the matrices Q, P1 (α ∗ )−1 Q, P2 (α ∗ )−1 Q and P3 (α ∗ )−1 Q for the fifth-order RadauIIA discretization of the FDE (4.1) (β = 1.8, s = 3, N = 28 − 1, τ = 0.1).
steps of the preconditioned GMRES method. The preconditioner P2 (α ∗ ) performs better than the preconditioners P1 (α ∗ ) and P3 (α ∗ ) in terms of CPU times. The numerical results show that the BTS-preconditioned GMRES performs better than the stationary BTS iteration. Furthermore, we note that the performance of the KPS-based polynomial preconditioners Pk (α ∗ ) (k = 1, 2, 3) is better than that of the BTS preconditioner tested in terms of computational times. In Fig. 4.1 we display the eigenvalues of the matrix Q before and after preconditioning with Pk (α ∗ ) (k = 1, 2, 3). These figures clearly show that the matrix without preconditioning is very ill conditioned while the matrices with the proposed preconditioning have tightly clustered eigenvalues and, thus, are well conditioned. In Fig. 4.2 we summarize results of the 2s − 1-order RadauIIA methods (RadauIIA-2s − 1) with s = 1, 2, 3, 4 for the FDE (4.1) with β = 1.5. Note that the first-order RadauIIA method is just the implicit Euler method. The spatial fractional operators of the FDE are discretized by WSGD scheme with mesh size h = 2−9 (N = 210 − 1). All the linear systems are solved by GMRES with preconditioner P2 (α ∗ ) (or Pˆ) except the ones associated with implicit Euler method, which are solved directly by the Gohberg-Semencul-type formula. We compare the accuracy against time step size and the total CPU time. Here, the error is measured in the discrete L2 -norm at t = 1. It is shown that higher-order RadauIIA methods require fewer steps to achieve a given accuracy. They are also faster methods provided that the given accuracy is high enough. From Fig. 4.2, we also see that the performance of the BTS preconditioner is better than the KPS-based preconditioner P2 (α ∗ ) in terms of computational cost for RadauIIA-3 method, while the preconditioner P2 (α ∗ ) performs better than the BTS preconditioner for RadauIIA-5 and RadauIIA-7 methods.
122
H. Chen, X. Wang and X. Li / Applied Mathematics and Computation 351 (2019) 116–123
Fig. 4.2. Accuracy versus time-step and CPU time for four RadauIIA methods applied to the FDE (4.1) with β = 1.5.
5. Conclusions In this paper we have proposed a KPS-based polynomial preconditioning strategy for systems which arise from implicit Runge-Kutta discretizations of the FDEs. The proposed preconditioner is constructed based on the KPS iteration method. One component of the KPS-based polynomial preconditioner has the same structure as the system derived from implicit Euler discretization of the FDEs. So we can reuse the high performance of implicit Euler discretization solver as the building block for our preconditioner. Theoretical and numerical results show that the proposed preconditioning method is effective and efficient. Acknowledgments The authors are very grateful to the referees for their constructive comments and valuable suggestions, which greatly improved the original manuscript of the paper. This work was supported by the National Natural Science Foundation of China (grant nos. 11301575 and 11471063), the Scientific and Technological Research Program of Chongqing Municipal Education Commission (grant nos. KJ1703055 and KJZD-M201800501), the Natural Science Foundation Project of CQ CSTC (nos. cstc2018jcyjAX0113 and cstc2018jcyjAX0794), and the Talent Project of Chongqing Normal University (grant no. 02030307– 0054). References [1] P. Amodio, L. Brugnano, A note on the efficient implementation of implicit methods for ODEs, J. Comput. Appl. Math. 87 (1997) 1–9. [2] O. Axelsson, R. Blaheta, R. Kohut, Preconditioning methods for high-order strongly stable time integration methods with an application for a DAE problem, Numer. Linear Algebra. Appl. 22 (2015) 930–949. [3] Z.Z. Bai, K.Y. Lu, J.Y. Pan, Diagonal and toeplitz splitting iteration methods for diagonal-plus-toeplitz linear systems from spatial fractional diffusion equations, Numer. Linear Algebra. Appl. (2017), doi:10.1002/nla.2093. [4] L. Brugnano, F. Iavernaro, C. Magherini, Efficient implementation of radau collocation methods, Appl. Numer. Math. 87 (2015) 100–113. [5] Z.H. Cao, Y. Wang, On validity of m-step multisplitting preconditioners for linear systems, Appl. Math. Comput. 126 (2002) 199–211. [6] H. Chen, A splitting preconditioner for the iterative solution of implicit Runge-Kutta and boundary value methods, BIT 54 (2014) 607–621. [7] H. Chen, Generalized kronecker product splitting iteration for the solution of implicit Runge-Kutta and boundary value methods, Numer. Linear Algebra. Appl. 22 (2015) 357–370. [8] H. Chen, Kronecker product splitting preconditioners for implicit runge-kutta discretizations of viscous wave equations, Appl. Math. Model. 40 (2016) 4429–4440. [9] H. Chen, A splitting preconditioner for implicit runge-kutta discretizations of a partial differential-algebraic equation, Numer. Algor. 73 (2016) 1037–1054. [10] H. Chen, W. Lv, T.T. Zhang, A kronecker product splitting preconditioner for two-dimensional space-fractional diffusion equations, J. Comput. Phys. 360 (2018) 1–14. [11] H. Chen, T.T. Zhang, W. Lv, Block preconditioning strategies for time-space fractional diffusion equations, Appl. Math. Comput. 337 (2018) 41–53. [12] G.J. Cooper, J.C. Butcher, An iteration scheme for implicit Runge-Kutta methods, IMA J. Numer. Anal. 3 (1983) 127–140. [13] M. Donatelli, M. Mazza, S. Serra-Capizzano, Spectral analysis and structure preserving preconditioners for fractional diffusion equations, J. Comput. Phys. 307 (2016) 262–279. [14] I. Gohberg, V. Olshevsky, Circulants, displacements and decompositions of matrices, Integral Equ. Oper. Theory 15 (1992) 730–743. [15] X.M. Gu, T.Z. Huang, H.B. Li, L. Li, W.H. Luo, On k-step CSCS-based polynomial preconditioners for toeplitz linear systems with application to fractional diffusion equations, Appl. Math. Lett. 42 (2015a) 53–58. [16] X.M. Gu, T.Z. Huang, X.L. Zhao, H.B. Li, L. Li, Strang-type preconditioners for solving fractional diffusion equations by boundary value methods, J. Comput. Appl. Math. 277 (2015b) 73–86. [17] E. Hairer, G. Wanner, Solving Ordinary Differential Equations II. Stiff and Differential Algebraic Problems, Springer-Verlag, Berlin, 1996. [18] Y.M. Huang, On m-step hermitian and skew-hermitian splitting preconditioning methods, J. Eng. Math. 93 (2015) 77–86. [19] L.O. Jay, T. Braconnier, A parallelizable preconditioner for the iterative solution of implicit Runge-Kutta type methods, J. Comput. Appl. Math. 111 (1999) 63–76. [20] S.L. Lei, H.W. Sun, A circulant preconditioner for fractional diffusion equations, J. Comput. Phys. 242 (2013) 715–725.
H. Chen, X. Wang and X. Li / Applied Mathematics and Computation 351 (2019) 116–123
123
[21] S.L. Lei, Y.C. Huang, Fast algorithms for high-order numerical methods for space-fractional diffusion equations, Int. J. Comput. Math. 94 (2017) 1062–1078. [22] X.L. Lin, M.K. Ng, H.W. Sun, A splitting preconditioner for toeplitz-like linear systems arising from fractional diffusion equations, SIAM J. Matrix Anal. Appl. 38 (2017) 1580–1614. [23] X.L. Lin, M.K. Ng, H.W. Sun, Efficient preconditioner of one-sided space fractional diffusion equation, BIT 58 (2018) 729–748. [24] K.A. Mardal, T.K. Nilssen, G.A. Staff, Order optimal preconditioners for implicit runge-kutta schemes applied to parabolic PDEs, SIAM J. Sci. Comput. 29 (2007) 361–375. [25] H. Moghaderi, M. Dehghan, M. Donatelli, M. Mazza, Spectral analysis and multigrid preconditioners for two-dimensional space-fractional diffusion equations, J. Comput. Phys. 350 (2017) 992–1011. [26] M. Ng, Iterative methods for Toeplitz systems, Oxford University Press, 2004. [27] I. Podlubny, Fractional differential equations, Math. Sci. Eng. 198 (1999). Academic Press, San Diego. [28] S. Perez-Rodriguez, S. Gonzalez-Pinto, B.P. Sommeijer, An iterated radau method for time-dependent PDEs, J. Comput. Appl. Math. 231 (2009) 49–66. [29] G.A. Staff, K.A. Mardal, T.K. Nilssen, Preconditioning of full implicit Runge-Kutta schemes for parabolic PDEs, Model. Identif. Control 27 (2006) 109–123. [30] W. Tian, H. Zhou, W. Deng, A class of second order difference approximation for solving space-fractional diffusion equations, Math. Comp. 84 (2015) 1703–1727. [31] P.J. van der Houwen, J.J.B. de Swart, Triangularly implicit iteration methods for ODE-IVP solvers, SIAM J. Sci. Comput. 18 (1997) 41–55. [32] J.V. Lent, S. Vandewalle, Multigrid methods for implicit runge-kutta and boundary value method discretizations of PDEs, SIAM J. Sci. Comput. 27 (2005) 67–92.