# A Banded Preconditioning Iteration Method for Time-Space Fractional Advection-Diffusion Equations

Authors: Min-Li Zeng; Guo-Feng Zhang
Journal: Mathematical Problems in Engineering
 (2017)

---

## Abstract

In this paper, we concentrate on the efficient solvers for the time-space fractional advection-diffusion equations. Firstly, the implicit finite difference schemes with the shifted Grünwald-Letnikov approximations for spatial fractional derivative and unshifted Grünwald-Letnikov approximations for time fractional derivative are employed to discretize time-space fractional advection-diffusion equations. The discretization results in a series of large dense linear systems. Then, a banded preconditioner is proposed and some theoretical properties for the preconditioning matrix are studied. Numerical implementations show that the banded preconditioner may lead to satisfactory experimental results when we choose appropriate bandwidth in the preconditioner.


---

## Body

## 1. Introduction

Consider the following initial-boundary value problem of time-space fractional advection-diffusion equations (TS-FADE):(1)D0Ctαux,t=d+x∂βux,t∂+xβ+d-x∂βux,t∂-xβ+v∂ux,t∂x+fx,t,xL≤x≤xR,t∈0,Tf,uxL,t=uxR,t=0,0≤t≤Tf,ux,0=0,x∈xL,xR,where f(x,t) is the source and sink term, the coefficient of anomalous diffusion satisfying d+(x)≥dmin+>0 and d-(x)≥dmin->0, the real constant v denotes the drift of the process, that is, the mean advection velocity, and β∈(1,2) and α∈(0,1) are the differentiation parameters.The Caputo time derivative of real orderα with α∈(0,1) in (1) is defined as [1, 2](2)D0Ctαux,t≔1Γ1-α∫0t∂ux,s/∂st-sαds.The right- and left-hand spatial fractional derivatives are defined by [3](3)∂βux,t∂+xβ=1Γ2-β∂2∂x2∫xLxuξ,tdξx-ξβ-1,∂βux,t∂-xβ=1Γ2-β∂2∂x2∫xxRuξ,tdξξ-xβ-1,respectively, where Γ(·) is the gamma function.The TS-FADE of form (1) arises in variety of research areas such as modeling chaotic dynamics of classical conservative systems [4], turbulent flow [5], groundwater contaminant transport [6], and applications in finance [7], image processing [8], physics [9], biological systems [10], and so on. Though many analytic approaches, such as the Fourier transform method, the Mellin transform method, and the Laplace transform methods, have been used to seek the closed-form solutions [11], there are few available analytical closed-form solutions for fractional differential equations (FDEs). Hence, it is important to find efficient methods to solve fractional differential equations.Traditional methods for solving FDEs generate the computational cost ofO(N3) and storage of O(N2) as a result of full coefficient matrices, where N is the grid point number. However, when we use the shifted Grünwald-Letnikov discretization scheme proposed by Meerschaert and Tadjeran [12] to approximate the FDE, we will obtain a special Toeplitz-like discretized coefficient matrix [13, 14]. Therefore, according to [13], the storage requirement is O(N) instead of O(N2), and the complexity of the matrix-vector multiplication only requires O(Nlog⁡N) operations by fast Fourier transform (FFT). By making use of this advantage, K. Wang and H. Wang [14] showed that the conjugate gradient normal residual (CGNR) method which has computational cost of O(Nlog2⁡N) is very efficient when the diffusion coefficients are very small; that is, the discretized systems are well-conditioned. However, when the diffusion coefficients are not small, the discretized diffusion coefficient matrices become ill-conditioned. Therefore, preconditioning technique must be used to improve the efficiency of the CGNR method. Many references proposed efficient methods. For example, Pang and Sun [15] proposed the multigrid (MG) method to solve the system. Lei and Sun [16] used the circulant preconditioned CGNR method, which is proved to be more efficient than the MG method. After that, Qu et al. [17] introduced a circulant and skew-circulant splitting method to solve the fractional diffusion equations. For more accuracy analysis of the discretization from the fractional difference equations, we refer to the documents [18–20]. Some efficient algorithms for the discretized linear systems from the space fractional difference equations can be found in [16, 21–24] and so on.In this paper, we will first discretize TS-FADE (1) by using the unshifted Grünwald-Letnikov approximation to the Caputo time fractional derivative and the shifted Grünwald-Letnikov approximation to the spatial fractional derivative. Then this discretization will lead to a system with Toeplitz-like coefficient matrix. Following the strategy proposed in [21], we will propose a banded preconditioning technique for solving the discretized system. Some theoretical results about the preconditioning matrix will be presented. Numerical experiments will also be given to testify the effectiveness of the new preconditioning technique.The remainder of this paper is organized as follows. In Section2, we give a detailed description for the discretization of TS-FADE (1). In Section 3, a banded preconditioner is constructed. Then some properties of the discretized system (16) are given and some spectral properties for the preconditioning matrix are studied. Numerical experiments are presented in Section 4 to illustrate the efficiency of the banded preconditioning technique for solving the TS-FADE. Finally, in Section 5, we draw some conclusions to end this paper.


## 2. Discretization of the TS-FADE

LetN and M be positive integers and h≔(xR-xL)/(N+1) and τ≔Tf/M be the sizes of spatial grid and time step, respectively. We first perform discretization in time using the unshifted Grünwald-Letnikov approximation to the Caputo time fractional derivative D0Ctαu(x,t). By making use of the abbreviations ui(m)≔u(xi,tm) and fi(m)≔f(xi,tm), where xi=xL+ih, tm=mτ,   i=0,1,…,N+1,  and m=0,1,…,M, we can obtain the following semidiscretized scheme:(4)TαM+1ux,t0ux,t1ux,t2⋮ux,tM=IM+1d+x∂βux,t0∂+xβ+d-x∂βux,t0∂-xβ+v∂ux,t0∂xd+x∂βux,t1∂+xβ+d-x∂βux,t1∂-xβ+v∂ux,t1∂xd+x∂βux,t2∂+xβ+d-x∂βux,t2∂-xβ+v∂ux,t2∂x⋮d+x∂βux,tM∂+xβ+d-x∂βux,tM∂-xβ+v∂ux,tM∂x+IM+1fx,t0fx,t1fx,t2⋮fx,tM,where IM+1∈R(M+1)×(M+1) is identity matrix and(5)
T
α
M
+
1
≔
τ
-
α
g
0
α
g
1
α
g
0
α
g
2
α
g
1
α
g
0
α
⋱
⋱
⋱
⋱
⋱
⋱
⋱
g
1
α
g
0
α
g
M
α
g
M
-
1
α
⋱
g
2
α
g
1
α
g
0
α
=
τ
-
α
g
0
α
0
⋮
T
α
M
τ
-
α
g
M
αrepresents the discrete Caputo derivative. Here gk(β) is defined as(6)g0β=1,gkβ=-1kk!ββ-1⋯β-k+1,fork=1,2,….Because of the zero initial condition, we can rewrite (4) as(7)TαMux,t1ux,t2⋮ux,tM=IMd+x∂βux,t1∂+xβ+d-x∂βux,t1∂-xβ+v∂ux,t1∂xd+x∂βux,t2∂+xβ+d-x∂βux,t2∂-xβ+v∂ux,t2∂x⋮d+x∂βux,tM∂+xβ+d-x∂βux,tM∂-xβ+v∂ux,tM∂x+IMfx,t1fx,t2⋮fx,tM,and here the first line of (7) has the form(8)g0αux,t1=d+x∂βux,t1∂+xβ+d-x∂βux,t1∂-xβ+v∂ux,t1∂x+fx,t1.We use the central difference approximation to ∂u(x,t1)/∂x and approximate the spatial derivative of order β(1<β≤2) by the shifted Grünwald-Letnikov approximations [12] for the right- and left-handed spatial fractional derivative; that is,(9)∂βuxi,tm∂+xβ=1hβ∑k=0i+1gkβui-k+1m+Oh,∂βuxi,tm∂-xβ=1hβ∑k=0N-i+2gkβui+k-1m+Oh.DenoteB≔(TαM)T,u(m)≔(u1m,u2m,…,uN(m))T,and f(m)≔(f1m,f2m,…,fNm)T,m=0,1,2,…,M.(10)D+≔diag⁡d+x1,d+x2,…,d+xN≡diag⁡d+,1,d+,2,…,d+,N,D-≔diag⁡d-x1,d-x2,…,d-xN≡diag⁡d-,1,d-,2,…,d-,N,(11)GβN≔-h-βg1βg0β0⋯00g2βg1βg0β0⋯0⋮g2βg1β⋱⋱⋮⋮⋱⋱⋱⋱0gN-1β⋱⋱⋱g1βg0βgNβgN-1β⋯⋯g2βg1β,(12)W≔0-10⋯010-1⋱⋮01⋱⋱0⋮⋱⋱⋱-10⋯010.According to the boundary condition u0(m)=uN+1(m)=0, the space-time discretization of (7) will lead to the following algebraic linear system in Kronecker form:(13)BT⊗IN+IM⊗Au=f~,where(14)A≔D+GβN+D-GβNT+v2hW,u≔[u(1);u(2);…;u(M)] using the MATLAB notation, and each u(m) is a vector of dimension N associated with the point m in time. Meanwhile, the right-hand side vector f~ in (13) is defined by f~≔[f(1);f(2);…;f(M)].Let(15)U≔u1,u2,…,uM,F≔f1,f2,…,fM.Then we can rewrite the linear system (13) into a matrix equation form(16)AU+UB=F,with A defined in (14) being a sum of some dense diagonal-times-Toeplitz matrices and B=(TαM)T being upper triangular matrix.Note thatA is a sum of some diagonal-times-Toeplitz matrices and B is upper triangular matrix, the right-hand side matrix in (16) is not low rank, and the size of B usually has much smaller size than that of A. Then the first column of U can be obtained by solving the shifted system (A+B11I)u(1)=f(1). All subsequent columns of U may be obtained with some substitutions as(17)A+BmmIum=fm-∑k=1m-1ukBkm,m=2,3,…,M,where B=(Bkm); that is, Bkm denotes the (k,m)th-entry of B.SinceBmm≡g0α=1,m=1,2,…,M. Then we have to solve these systems with the same coefficient matrix A+I=D+GβN+D-(GβN)T+(v/2h)W+I, where D+ and D- are positive diagonal matrices, I is an identity matrix, and W and GβN are described previously.


## 3. Banded Preconditioning Iteration Method

Define(18)A≔BmmI+A=I+D+GβN+D-GβNT+ηW,where η=v/2h<0 and the matrices D+, D-, W, and GβN are defined in Section 2. It can be easily found that W is a skew-symmetric matrix.Obviously, the main work for solving the discretized system (17) concentrates on solving the system with coefficient matrix A for m=1,2,3,…,M. Since A is a nonsymmetric matrix, then we can use the Krylov subspace methods [25], such as the generalized minimal residual (GMRES) method, to solve the linear systems (17). The preconditioning techniques are often employed to improve the performance and reliability of the Krylov subspace methods.Denote(19)Gβ,kN≔-h-βg1βg0β⋮g1βg0βgkβ⋱⋱⋱⋱⋱gkβ⋯g1β,2≤k≤N,and define a banded preconditioner Pk as(20)Pk≔I+D+Gβ,kN+D-Gβ,kNT+ηW.It can be easily seen that Gβ,kN is a banded Toeplitz matrix with bandwidth k+1 containing the central diagonals of GβN. Therefore, Pk is a banded matrix with bandwidth 2k-1.To study the property of the preconditionerPk, we will give the following lemmas first.Lemma 1 (see [12]).
Let1<β≤2 and gj(β) be defined in (6). Then we have(21)g0β=1,g1β=-β<0,g2β>g3β>⋯>0,∑j=0∞gjβ=0,∑j=0ngjβ<0,forn≥1,gjβ=Oj-β+1.Lemma 2.
For1<β≤2, the matrix Gβ,kN defined in (19) is a strongly diagonally dominant matrix and the real part of the spectrum of Gβ,kN is contained in the open interval (0,2h-ββ).Proof.
According to Lemma1, it holds that ∑j=0ngj(β)<0forn≥1. Therefore, for each row of the matrix Gβ,kN, we have ∑j=0lgj(β)<0 or, equivalently, g1β=-g1(β)>∑j=0,j≠1lgj(β) for l=1,2,…,k. That is, g1β≥∑j=0,j≠1lgjβ(l=1,2,…,k). Hence, Gβ,kN is a strongly diagonally dominant matrix.
Moreover, by simple computation, we can obtain the first column of the matrixHk≔(1/2)(Gβ,kN+(Gβ,kN)T) as(22)-h-βg1β,12g0β+g2β,12g3β,…,12gkβ,0,…,0T︸N.Hence, the Gershgorin circles of Hk are centered at -h-βg1(β)=βh-β. Moreover, according to Lemma 1, the largest radius rmax is at most(23)rmax=h-βg0β+g2β+g3β+⋯+gkβ=h-β∑j=0kgjβ-g1β<h-βg1β=h-ββ.It then follows that the real part of the spectrum of Gβ,kN is contained in the open interval (0,2h-ββ).Next, we will concentrate on the uniqueness of the solution of the matrix equation (16).Theorem 3.
LetG≔hβ·(GβN+(GβN)T), where GβN is defined as in (11). Then G is a symmetric positive definite matrix.Proof.
Suppose that{λj(G),j=1,2,…,N} are the eigenvalues set of the matrix G. According to the Gershgorin theorem, all eigenvalues of G are contained within a complex disk centered at -2g1(β)=2β, with the maximum radius being in line ⌊N/2⌋+1 or in line ⌊N/2⌋ (dependent on whether N is odd or even). Therefore, the maximum radius rmax satisfies(24)rmax=2∑k=0,k≠1N/2+1gkβ<-2g1β=2β.That is, {Re⁡(λj(G))}⊂(0,4β).
BecauseG is a real symmetric matrix, then the eigenvalues are all real. So G is symmetric positive definite.Theorem 4.
Suppose thatF defined in (16) satisfies F≠0, and then there exists only one solution for the Sylvester equation (16).Proof.
From [26, 27], we know that when C≠0, the continuous Sylvester equation (16) has a unique solution if there is no common eigenvalue between A and -B. It is easy to know that B is upper triangular matrix and its eigenvalues are all positive. Note that A+AT=(D+GβN+(D+GβN)T)+(D-(GβN)T+(D-GβN)T)T+(v/2h)W+WT=D+GβN+GβN)T+D-(GβN)T+GβN)≔D1+D2. According to Theorem 3, if the real parts of the eigenvalues of GβN are positive, then the real parts of all eigenvalues of D+GβN are also positive. Furthermore, the matrices D1 and D2 are both symmetric positive definite. Therefore, A+AT is symmetric positive definite, and the real parts of the eigenvalues of A are positive. Therefore, there is no common eigenvalues between the matrices A and -B.Theorem 5.
For2≤k≤N, the preconditioner Pk defined in (20) is nonsingular.Proof.
SinceD+ is a positive diagonal matrix and from Lemma 2, the real part of the spectrum of Gβ,kN is positive. Then the real part of the spectrum of D+Gβ,kN is also positive. Using the same strategy, we can obtain that the real part of the spectrum of D-(Gβ,kN)T is also positive.
Assume thatPk is a singular matrix and there exists a nonzero eigenvector x such that Pkx=0. Then it follows that(25)I+D+Gβ,kN+D-Gβ,kNTx=-ηWx.Multiplying x∗/x∗x by the left of the above equation, we get(26)x∗I+D+Gβ,kN+D-Gβ,kNTxx∗x=-ηx∗Wxx∗x.As W is a skew-symmetric matrix and η<0 is a real number, then -η(x∗Wx/x∗x) is a pure imaginary. Because the real part of x∗(D+Gβ,kN+D-(Gβ,kN)T)x/x∗x is a positive number, then the real part of x∗(I+D+Gβ,kN+D-(Gβ,kN)T)x/x∗x is a positive number. According to (26), we can easily find that there is a contradiction because the real part of the right hand side of (26) is zero. That is, if Pkx=0, it must hold that x=0; that is, the matrix Pk is nonsingular.Theorem 6.
For2≤k≤N and 1<β≤2, the relative difference between A and Pk can be described as(27)A-Pk∞A∞≤Ok-β,which implies Pk can be an efficient preconditioner for the coefficient matrix A as k increases.Proof.
The proof process is similar to [13].At the end of this section, we will draw the main steps for solving the discretized system (13).Step 1.
Compute the matricesA, B, and F in (16).Step 2.
Solve the systemg0(α)u(1)+Au(1)=f(1) by using the banded preconditioner Pk.Step 3.
Computeb(m)=f(m)-∑k=1m-1u(k)Bkm.Step 4.
Form=2,3,…,M, solve the systems (A+BmmI)u(m)=b(m) by using the banded preconditioner Pk.Step 5.
Obtain the numerical solutionU≔[u(1),u(2),…,u(M)].


## 4. Numerical Results

In this section, we will test the feasibility and effectiveness of the banded preconditioning iteration method with the preconditionerPk for TS-FADE (1). All runs at each time step are started from zero initial guess and terminated if the current iteration satisfies ERR≤10-7 or the elapsed CPU time is exceeding 2000 in seconds, where(28)ERR≔bm-Aum,k2bm-Aum,02.All the experiments are performed in MATLAB R2012a on Intel® Core™ i7-3770 CPU 3.40 GHz and 8.00 GB of RAM, with machine precision 10-16. The average iteration counts (denoted by “Iter”), elapsed CPU time in seconds (denoted by “CPU”), and total iteration counts for solving the discretized system (17) (denoted by “IT”) are reported. Here “Iter” is defined as(29)Iter=ITM.Example 7.
Consider TS-FADE (1) with xL=0, xR=2, and Tf=1. The coefficient functions are(30)d+x=1+x4,d-x=1+2-x4,and the forcing function is(31)fx,t=4π2αsin⁡t2π-4sin⁡t2πΓ5Γ5-βd+xx4-β+d-x2-x4-β-4Γ4Γ4-βd+xx3-β+d-x2-x3-β+4Γ3Γ3-βd+xx2-β+d-x2-x2-β-16sin⁡t2πx3-3x2+2x.Then the true solution for the corresponding TS-FADE is u(x,t)=4sin⁡((t/2)π)x2(2-x)2.This example is a modification of the example from [13, 15, 23]. By applying the unshifted Grünwald-Letnikov approximations to discretize the fractional time derivative, the shifted Grünwald-Letnikov approximations to the fractional spatial derivative, and the central difference finite scheme to discretize the first-order derivative of spatial derivative, we can obtain the linear systems (17) step by step. To obtain the numerical solutions of this problem, we have to solve linear equations with the same coefficient matrix and different right-hand sides in each time step. In our experiments, we choose the parameters η=-1. To be more persuasive, we test the methods according to different choices of the derivative parameters α∈{0.4,0.8} andβ∈{1.3,1.6,1.9} and different mesh grids N+1∈{27,28,29,210,211,212,213} andM=24.When we use preconditioned GMRES method, we test the banded preconditioners with varying bandwidth (denoted as “P(k)-GMRES”), the Strang-like circulant preconditioner (denoted as “S-GMRES”), the R. Chan-like circulant preconditioner (denoted as “R-GMRES”), and the T. Chan-like circulant preconditioner [28] (denoted as “T-GMRES”):(32)I+d¯+C+d¯-CT+ηWc,where d¯+ and d¯- are the mean values of the diagonals of the diagonal matrices D+ and D-, respectively. The matrix Wc is defined by (33)0-10⋯110-1⋱⋮01⋱⋱0⋮⋱⋱⋱-1-1⋯010.The matrix C is chosen as Strang’s circulant matrix (Cs), R. Chan’s circulant matrix (Cr), and T. Chan’s circulant matrix (Ct), respectively. Further, the entries of the first column of Strang’s circulant matrix Cs are given by(34)-h-βg1β,g2β,…,gN/2β,0,…,0,g0βT.The entries of the first column of R. Chan’s circulant matrix Cr are given by(35)-h-βg1β,g2β,…,gN-1β,gNβ+g0βT.The entries of the first column of T. Chan’s circulant matrix Ct are given by(36)-h-βg1β,N-1g2β+gNβN,…,2gN-1β+N-2g2βN,gNβ+N-1g0βNT.We test the preconditioners P(4) and P(16), which means the half width of the corresponding banded preconditioner is (N+1)/4 and (N+1)/16, respectively.We show the experimental results forα=0.4 and α=0.8 in Tables 1 and 2, respectively. When the elapsed CPU time is exceeding 2000 in seconds, the results will be denoted by “—”; from these tables, we see that the performance for the banded preconditioner with the half width being 1/16 of the problem size is better compared to the circulant preconditioners. However, if the half width is 1/4 of the problem size, the performance reaches the best.Table 1
Numerical results of the tested methods forα=0.4 and M=24.
β N + 1 S-GMRES T-GMRES R-GMRES P(16)-GMRES P(4)-GMRES IT CPU IT CPU IT CPU IT CPU IT CPU 1.3 2 7 18 0.18 22 0.24 18 0.15 8 0.07 5 0.05 2 8 20 0.51 26 0.58 19 0.43 8 0.18 5 0.12 2 9 21 1.60 30 2.23 21 1.61 8 0.63 5 0.53 2 10 22 9.12 35 14.14 22 9.14 8 3.02 5 2.33 2 11 23 51.10 39 83.54 23 51.21 8 18.57 5 12.44 2 12 23 306.29 44 560.60 23 307.80 8 110.39 5 75.95 2 13 — — — — — — 8 741.27 5 495.16  1.6 2 7 11 0.09 18 0.17 11 0.11 8 0.08 4 0.06 2 8 12 0.28 22 0.46 12 0.28 8 0.17 4 0.11 2 9 12 0.97 26 1.93 12 0.97 8 0.63 4 0.43 2 10 12 5.26 30 12.11 12 5.30 8 3.46 4 2.32 2 11 12 28.88 36 77.01 12 28.60 8 18.7 4 12.4 2 12 12 169.48 43 547.79 12 169.37 8 109.46 4 73.25 2 13 12 1153.13 — — 12 1157.55 8 742.61 4 496.21  1.9 2 7 8 0.08 18 0.15 8 0.08 6 0.06 3 0.04 2 8 8 0.19 22 0.46 8 0.19 6 0.14 3 0.1 2 9 8 0.71 27 2.00 8 0.70 6 0.49 3 0.36 2 10 8 3.79 34 13.58 8 3.81 6 2.71 3 1.93 2 11 8 21.18 44 93.08 8 20.57 6 14.45 3 10.36 2 12 8 121.58 57 717.37 8 121.41 6 85.41 3 61.08 2 13 7 745.34 — — 7 742.52 6 578.17 3 413.03Table 2
Numerical results of the tested methods forα=0.8 and M=24.
β N + 1 S-GMRES T-GMRES R-GMRES P(16)-GMRES P(4)-GMRES IT CPU IT CPU IT CPU IT CPU IT CPU 1.3 2 7 18 0.16 22 0.20 18 0.16 8 0.07 5 0.05 2 8 20 0.43 26 0.55 19 0.40 8 0.17 5 0.12 2 9 21 1.59 30 2.22 21 1.60 8 0.63 5 0.43 2 10 22 9.12 35 14.03 22 9.11 8 3.04 5 2.32 2 11 23 51.60 39 83.07 23 51.10 8 18.52 5 12.83 2 12 23 302.65 44 559.36 23 308.72 8 109.32 5 73.38 2 13 — — — — — — 8 736.5 5 493  1.6 2 7 11 0.10 18 0.15 11 0.10 8 0.07 4 0.04 2 8 12 0.27 22 0.46 12 0.37 8 0.17 4 0.27 2 9 12 0.96 26 1.92 12 0.97 8 0.62 4 0.44 2 10 12 5.30 30 12.02 12 5.30 8 3.42 4 2.32 2 11 12 28.66 36 77.11 12 30.14 8 18.55 4 12.38 2 12 12 169.48 43 547.61 12 170.06 8 109.37 4 73.19 2 13 12 1146.07 — — 12 1146.50 8 738.33 4 494.49  1.9 2 7 8 0.07 18 0.15 8 0.07 6 0.06 3 0.04 2 8 8 0.19 22 0.46 8 0.19 6 0.14 3 0.1 2 9 8 0.70 27 2.00 8 0.72 6 0.49 3 0.35 2 10 8 4.22 34 13.76 8 3.82 6 2.7 3 1.94 2 11 8 20.79 44 93.85 8 20.52 6 14.44 3 10.39 2 12 8 121.68 57 717.90 8 121.76 6 85.35 3 61.15 2 13 7 740.29 — — 7 742.03 6 576.51 3 411.58


## 5. Conclusions

In this paper, we concentrate on the efficient solvers for time-space fractional advection-diffusion equations. By employing the implicit finite difference schemes with the shifted Grünwald-Letnikov approximations for spatial fractional derivative and unshifted Grünwald-Letnikov approximations for time fractional derivative, we can discretize the time-space fractional advection-diffusion equations. To solve the discretization, we have to solve a series of linear systems with the same coefficient matrix. Therefore, we proposed a banded preconditioning iteration method for the numerical solution of the resulting linear systems. Theoretical properties for the preconditioning matrix are studied in detail. Numerical implementations show that the banded preconditioner leads to satisfactory experimental results when we choose appropriate bandwidth in the preconditioner.



---

*Source: 6961296-2017-01-24.xml*