A Domain Decomposition Method for Time Fractional Reaction-Diffusion Equation

Abstract

The computational complexity of one-dimensional time fractional reaction-diffusion equation is 𝑂(𝑁²𝑀) compared with 𝑂(𝑁𝑀) for classical integer reaction-diffusion equation. Parallel computing is used to overcome this challenge. Domain decomposition method (DDM) embodies large potential for parallelization of the numerical solution for fractional equations and serves as a basis for distributed, parallel computations. A domain decomposition algorithm for time fractional reaction-diffusion equation with implicit finite difference method is proposed. The domain decomposition algorithm keeps the same parallelism but needs much fewer iterations, compared with Jacobi iteration in each time step. Numerical experiments are used to verify the efficiency of the obtained algorithm.

1. Introduction

Fractional equations can be used to describe some physical phenomenon more accurately than the classical integer order differential equation. The reaction-diffusion equations play an important role in dynamical systems of mathematics, physics, chemistry, bioinformatics, finance, and other research areas. There has been a wide variety of analytical and numerical methods proposed for fractional equations [1–7], for example, finite difference method [8], finite element method [9], Adomian decomposition method [10], and spectral technique [11]. Interest in fractional reaction-diffusion equations has increased [12].

Domain decomposition methods (DDM) solve a boundary value problem by splitting it into smaller boundary value problems on subdomains and iterating it to coordinate the solution between adjacent subdomains [13]. A coarse problem with one or few unknowns per subdomain is used to further coordinate the solution between the subdomains globally. The DDM can be divided into two categories: the overlapping and nonoverlapping [14]. Chan and Mathew [15] gave a survey on iterative domain decomposition techniques that had been developed for solving several kinds of partial differential equations, including elliptic, parabolic, and differential systems such as the Stokes problem and mixed formulations of elliptic problems. The problems on the subdomains are almost independent, which makes domain decomposition methods suitable for parallel computing. Parallel computing is used to solve intensive computation applications simultaneously [16], such as particle transport [17, 18] and fast multipole methods [19]. It is time consuming to numerically solve fractional differential equations for long time tail. Parallel computing [20–22] can be used to overcome the computational challenge of fractional approximation. DDM will embody large potential for a parallelization of the numerical solution for fractional equations. Until today the power of DDM for approximating fractional derivatives and solving fractional differential equations has not been recognized.

This paper focuses on the Caputo fractional reaction-diffusion equation: ^𝐶₀𝐷^𝛼_𝑡𝑢()()=𝜕𝑥,𝑡+𝜇𝑢𝑥,𝑡₂𝑢(𝑥,𝑡)𝜕𝑥²(),(),𝑢()(𝑥)+𝐾𝑓𝑥,𝑡0<𝛼<1𝑥,0=𝑔,𝑥∈0,𝑥_𝑅,𝑢(0,𝑡)=𝑢𝑥_𝑅,𝑡=0,𝑡∈[0,𝑇] (1) on a finite domain 0≤𝑥≤𝑥_𝑅 and 0≤𝑡≤𝑇. The 𝜇>0 and 𝐾 are constants. If 𝛼 equals 1, (1) is the classical reaction-diffusion equation. The fractional derivative is in the Caputo form.

2. Background

2.1. Numerical Solution

The fractional derivative of 𝑓(𝑡) in the Caputo sense is defined as [23] ^𝐶₀𝐷^𝛼_𝑡𝑓(𝑡)=1Γ(1−𝛼)^𝑡₀𝑓^′(𝜉)(𝑡−𝜉)^𝛼𝑑𝜉,(0<𝛼<1). (2)

If 𝑓^′(𝑡) is continuous bounded derivatives in [0,𝑇] for every 𝑇>𝑎, we can get ₀𝐷_𝑡^{𝛼𝑓(𝑡)}=lim_{𝜉→0,𝑛𝜉=𝑡}𝜉^𝛼𝑛_𝑖=0(−1)^𝑖𝛼𝑖=𝑓(0)𝑡^−𝛼Γ()+11−𝛼Γ(1−𝛼)^𝑡₀𝑓^′(𝜉)(𝑡−𝜉)^𝛼𝑑𝜉. (3)

Define 𝜏=𝑇/𝑁, ℎ=𝑥_𝑅/(𝑀+1), 𝑡_𝑛=𝑛𝜏, and 𝑥_𝑖=0+𝑖ℎ for 0≤𝑛≤𝑁, 0≤𝑖≤𝑀+1. Define 𝑢^𝑛_𝑖, 𝑓^𝑛_𝑖, and 𝑔_𝑖 as the numerical approximation to 𝑢(𝑥_𝑖,𝑡_𝑛), 𝑓(𝑥_i,𝑡_𝑛), and 𝑔(𝑥_𝑖). We can get [12] ^𝐶₀𝐷^𝛼_𝑡𝑢(𝑥,𝑡)|^𝑡𝑛_𝑥𝑖=1()×𝑏𝜏Γ1−𝛼₀𝑢^𝑛_𝑖−^𝑛−1_𝑘=1𝑏_{𝑛−𝑘−1}−𝑏_𝑛−𝑘𝑢^𝑘_𝑖−𝑏_𝑛−1𝑢⁰_𝑖+○𝜏^2−𝛼, (4) where 1≤𝑖≤𝑀, 𝑛≥1, and 𝑏_𝑙=𝜏^1−𝛼)1−𝛼(𝑙+1^1−𝛼−𝑙^1−𝛼,𝑙≥0. (5)

By using center difference scheme for 𝜕²𝑢(𝑥,𝑡)/𝜕𝑥², we can get 𝜕²𝑢(𝑥,𝑡)𝜕𝑥^2𝑡𝑛_𝑥𝑖=1ℎ²𝑢^𝑛_𝑖+1−2𝑢^𝑛_𝑖+𝑢^𝑛_𝑖−1+○ℎ². (6)

The implicit finite difference approximation for (1) is 1()𝑏𝜏Γ1−𝛼⁰𝑢^𝑛_𝑖−^𝑛−1_𝑘=1𝑏_{𝑛−𝑘−1}−𝑏_𝑛−𝑘𝑢^𝑘_𝑖−𝑏_𝑛−1𝑢⁰_𝑖+𝜇𝑢^𝑛_𝑖=𝑢^𝑛_𝑖+1−2𝑢^𝑛_𝑖+𝑢^𝑛_𝑖−1ℎ²+𝐾𝑓^𝑛_𝑖. (7)

Define 𝑠=2/ℎ²+𝑏₀𝜏⁻¹/Γ(1−𝛼)+𝜇, 𝑈^𝑛=(𝑢^𝑛₁,𝑢^𝑛₂,…,𝑢^𝑛_𝑀)^𝑇, 𝐹^𝑛=(𝑓^𝑛₁,𝑓^𝑛₂,…,𝑓^𝑛_𝑀)^𝑇, and 𝑟_𝑙 as 𝑟_𝑙=𝑏_𝑙−𝑏_𝑙+1𝑠. (8)

Equation (7) evolves as 𝐴𝑈^𝑛=^𝑛−1_𝑘=1𝑟_{𝑛−1−𝑘}𝑈^𝑘+𝑏_𝑛−1𝑈⁰+𝐾𝐹^𝑛, (9) where matrix 𝐴 is a tridiagonal matrix, defined by 𝐴_𝑀×𝑀=1𝑠−ℎ²−1ℎ²1𝑠−ℎ²1⋅⋅⋅⋅⋅−ℎ²−1ℎ²𝑠 (10)

Because 𝜇>0, 𝑏₀>0, the elements of matrix 𝐴 satisfy |𝑠|>|−1/ℎ²|+|−1/ℎ²|. This means that matrix 𝐴 is strictly diagonally dominant.

2.2. Computational Challenge

In order to get 𝑈^𝑛, the right-sided computation of (9) should be performed and tridiagonal linear system should be solved. There are mainly many constant vector multiplications and many vector vector additions in the right-sided computation.(1)The constant vector multiplications are 𝑉^′=𝑏_𝑛−1𝑈⁰, 𝑉^𝑘=𝑟_{𝑛−1−𝑘}𝑈^𝑘, and 𝑉^′′=𝐾𝐹^𝑛.(2)The vector vector additions are 𝑉=𝑉^′+∑^𝑛−1_𝑘=1𝑉^𝑘+𝑉^′′.(3)After solving tridiagonal linear system 𝐴𝑈^𝑛+1=𝑉, we get 𝑈^𝑛+1.

The Thomas algorithm for tridiagonal systems needs 5𝑀 multiplications and 3𝑀 additions. The computational complexity of 𝐴𝑈^𝑛=𝑉 is 𝑂(𝑀). The total computation of (9) is determined by ∑^𝑛−1_𝑘=1𝑟_{𝑛−1−𝑘}𝑈^𝑘, which means (𝑛−1)𝑀 multiplications and (𝑛−2)𝑀 additions for each time step; ^𝑁_𝑛=1(2𝑛𝑀−3𝑀)=𝑂𝑁²𝑀. (11) The computational complexity of (1) is 𝑂(𝑁²𝑀), while the computational complexity of classical one-dimensional reaction-diffusion equation is only 𝑂(𝑁𝑀). The computational cost of (11) varies linearly along the number of grid points but squares with the number of time steps.

3. Domain Decomposition Method

3.1. DDM with Two Subdomains

Similar to the classical alternating Schwarz method [13, 24], the domain Ω=[0,𝑥_𝑅]=𝑝₀,𝑝₁,…,𝑝_𝑀 can be divided into two subdomains Ω_𝑎 and Ω_𝑏. There are 𝑀+1 grid points for Ω. Ω_𝑎={𝑝₀,𝑝₁,…,𝑝_𝑚,𝑝_𝑚+1} and Ω_𝑏={𝑝_𝑚,𝑝_𝑚+1,…,𝑝_𝑀}, where 0<𝑚<𝑀. The global physical boundary is defined in (1). 𝜕Ω_𝑎 (the right boundary of Ω_𝑎) and 𝜕Ω_𝑏 (the left boundary of Ω_𝑏) are called artificial internal boundary.

In order to approximate the time fractional equation on the two subdomains separately, the following iterative procedure can be performed. For each time step, the right hand side of (9) is calculated at first and the 𝑈^𝑛_Ω𝑖 is given as initial guess 𝑈_Ω^𝑛−1_𝑖, where 𝑖=𝑎,𝑏. The better approximation of 𝑈^𝑛_Ω𝑖 can be obtained iteratively. During each iteration, which is inside of a time step, the time fractional equation is solved in the subdomain Ω_𝑎, using the approximation of the previous iteration from Ω_𝑏 on 𝜕Ω_𝑎 as follows: ^𝐶₀𝐷^𝛼_𝑡𝑢()()=𝜕𝑥,𝑡+𝜇𝑢𝑥,𝑡₂𝑢(𝑥,𝑡)𝜕𝑥²(),(),𝑢()(𝑥)+𝐾𝑓𝑥,𝑡0<𝛼<1𝑥,0=𝑔,𝑥∈0,𝑥_𝑅,𝑢(𝑥,0)=𝑢^𝑛_{𝑏,𝑚+1,previous},𝑥on𝜕Ω_𝑎,𝑢0,𝑡_𝑛=0, (12) where 𝑢^𝑛_{𝑏,𝑚+1,previous} stands for the previous solution of grid point 𝑝_𝑚+1 in subdomain Ω_𝑏. The better approximation 𝑈^𝑛_𝑎,new is obtained. 𝑈^𝑛_𝑎,new is defined as {𝑢^𝑛_𝑎,1,new,…,𝑢^𝑛_{𝑎,𝑚,new}}. 𝑈^𝑛_{𝑎,previous} is defined as {𝑢^𝑛_{𝑎,1,previous},…,𝑢^𝑛_{𝑎,𝑚,previous}}. The definitions for 𝑈^𝑛_{𝑏,previous} and 𝑈^𝑛_𝑎,new are similar.

Then, we solve the time fractional equation within the subdomain of Ω_𝑏, using the approximation of the previous iteration from Ω_𝑎 on 𝜕Ω_𝑏 as follows: ^𝐶₀𝐷^𝛼_𝑡𝑢()()=𝜕𝑥,𝑡+𝜇𝑢𝑥,𝑡₂𝑢(𝑥,𝑡)𝜕𝑥²(),(),𝑢()(𝑥)+𝐾𝑓𝑥,𝑡0<𝛼<1𝑥,0=𝑔,𝑥=𝑥_𝑅𝑢(𝑥,0)=𝑢^𝑛_{𝑎,𝑚,previous},𝑥on𝜕Ω_𝑏𝑢𝑥_𝑅,𝑡_𝑛=0, (13) where 𝑢^𝑛_{𝑎,𝑚,previous} stands for the previous solution of grid point 𝑝_𝑚 in subdomain Ω_𝑎.

The two local time fractional equations in Ω_𝑎 and Ω_𝑏 are connected by the artificial boundary condition. The artificial boundary condition on the internal boundary 𝜕Ω_𝑎 of subdomain Ω_𝑎 is provided by 𝑢^𝑛_{𝑏,𝑚+1,previous} from subdomain Ω_𝑏, and vice versa. The approximation 𝑢^𝑛_{𝑎,𝑚,previous} and 𝑢^𝑛_{𝑏,𝑚+1,previous} may change until converged to the true solution. So, in an inner iteration of each time step, the two time fractional equations need to exchange two sets of data (send one and receive one) to update the artificial boundary conditions.

3.2. A Domain Decomposition Algorithm

Section 3.1 shows the procedure of DDM for time fractional equation with two subdomains. It is not hard to extend the method of Section 3.1 to more than two subdomains. The domain Ω can be decomposed into a set of 𝑃 subdomains {Ω_𝑝}^𝑃_𝑝=1 with Ω=∪^𝑃_𝑝=1Ω_𝑝. For time step 𝑛, Ω₁ has one global boundary 𝑥=0 and one artificial inner boundary 𝜕Ω_1,𝑏. Ω_𝑃 has one global boundary 𝑥=𝑥_𝑅 and one artificial inner boundary 𝜕Ω_𝑃,𝑎. The Ω_𝑝 (1<𝑝<𝑃) has two artificial inner boundaries 𝜕Ω_𝑃,𝑎 and 𝜕Ω_𝑝,𝑏. Ω_𝑝∩Ω_𝑝+1≠Φ means that the neighboring subdomains have explicit overlap.

The iterative procedure for the time step 𝑛+1 is similar to Section 3.1. The current iteration of Ω_𝑝 uses the data of previous iteration of its neighboring subdomains. Assuming 𝑀 is divisible with 𝑃, the domain decomposition algorithm is shown in Algorithm 1.

(1) [ℎ!]   input: 𝑀,𝑃,𝑥𝑅,𝑇,𝜇,𝐾,𝜖 et al.  Output: 𝑈 (2) 𝐼←𝑀/𝑃 (3) Allocate memory space 𝑈𝑁𝐼,𝑃,𝐴𝐼,3,𝑃,𝐹𝑁𝐼,𝑃,𝑉3𝐼,𝑃 et al. (4) Init matrices 𝑈,𝐴,𝐹,𝑏,𝑟 et al. (5) Declare local variables 𝛿𝑃,𝛿 (6) 𝑉1→31→𝐼,1→𝑃←0 (7) Get 𝑈01→𝐼,1→𝑃 with initial boundary. (8) 𝑡𝑜𝑡𝑎𝑙𝐼𝑡𝑒𝑟𝑎𝑡𝑖𝑜𝑛←0 (9) for  𝑛=0  to  𝑁−1 step by 1 do (10)  for 𝑝=1  to  𝑃 step by 1 do (11)     𝑈𝑛+11→𝐼,𝑝←𝑏𝑛𝑈01→𝐼,𝑝+𝐾𝐹𝑛+11→𝐼,𝑝 (12)     for  𝑘=1  to  𝑛 step by 1 do (13)     𝑈𝑛+11→𝐼,𝑝←𝑈𝑛1→𝐼,𝑝+𝑈𝑘1→𝐼,𝑝𝑟𝑛−𝑘 (14)     𝛿←1.0 (15)     𝑙𝑜𝑐𝑎𝑙𝐼𝑡𝑒𝑟𝑎𝑡𝑖𝑜𝑛←0 (16)     while 𝛿>𝜖 do (17)     for 𝑝=1  to  𝑃 step by 1 do (18)        𝑉21→𝐼,𝑝←𝑈𝑛+11→𝐼,𝑝 (19)        if 𝑝>1 then (20)        𝑉21,𝑝←𝑉21,𝑝+𝑉1𝐼,𝑝−1/ℎ2 (21)        if 𝑝<𝑃 then (22)        𝑉2𝐼,𝑝←𝑉2𝐼,𝑝+𝑉1𝐼,𝑝+1/ℎ2 (23)     for  𝑝=1  to  𝑃 step by 1 do (24)        solve 𝐴1→𝐼,3,𝑝𝑉31→𝐼,𝑝=𝑉21→𝐼,𝑝 (25)        𝛿𝑝=max𝑉1𝑖,𝑝−𝑉3𝑖,𝑝𝐼𝑖=1 (26)        𝛿=max𝛿,𝛿𝑝 (27)     𝑙𝑜𝑐𝑎𝑙𝐼𝑡𝑒𝑟𝑎𝑡𝑖𝑜𝑛←𝑙𝑜𝑐𝑎𝑙𝐼𝑡𝑒𝑟𝑎𝑡𝑖𝑜𝑛+1 (28)     𝑡𝑜𝑡𝑎𝑙𝐼𝑡𝑒𝑟𝑎𝑡𝑖𝑜𝑛←𝑡𝑜𝑡𝑎𝑙𝐼𝑡𝑒𝑟𝑎𝑡𝑖𝑜𝑛+𝑐𝑜𝑢𝑛𝑡 (29)     for 𝑝=1  to 𝑃 step by 1 do (30)     𝑈𝑛+11→𝐼,𝑝←𝑉31→𝐼,𝑝 (31) Output the information (32) Free memory space

Algorithm 1 

Domain decomposition algorithm for time fractional reaction-diffusion equation.

In Algorithm 1, there are some fast algorithms to solve the tridiagonal matrix 𝐴_{1→𝐼,3,𝑝}𝑉³_{1→𝐼,𝑝}=𝑉²_{1→𝐼,𝑝}, such as Thomas algorithm. 𝜖 is a threshold, such as 10⁻⁶. The signal 𝑙𝑜𝑐𝑎𝑙𝐼𝑡𝑒𝑟𝑎𝑡𝑖𝑜𝑛 is used to count how many iterations are needed in each time step. The data exchange between neighboring iterations is shown in lines 19–22. From the view of computer science, lines 2–6, lines 7–30, and lines 31-32 are preprocessing procedure, numerical solver, and postprocessing procedure.

3.3. Analysis

The presented DD algorithm updates the artificial boundary condition in a Jacobi fashion, using approximation from all the relevant neighboring subdomains from the previous iteration for each time step. A subdomain only exchanges two sets of data for one artificial boundary with its neighbor. Therefore, the subdomain solved in Algorithm 1 can be carried out almost completely independently, thus making the method inherently as parallel as the Jacobi iteration. The DD algorithm keeps the good parallelism of Jacobi iteration but needs fewer inner iterations in each time step; see Section 4. Equation (9) can be regarded as approximation of a special integer order reaction-diffusion equation. The stability and convergence analysis of integer order reaction-diffusion equation can refer to Mathew’s book [25].

4. Numerical Example

The following Caputo fractional reaction-diffusion equation [12] was considered, as shown in (14): ^𝐶₀𝐷^𝛼_𝑡𝑢()()=𝜕𝑥,𝑡+𝜇𝑢𝑥,𝑡₂𝑢(𝑥,𝑡)𝜕𝑥²(),(),𝑢()(),𝑢()()+𝐾𝑓𝑥,𝑡0<𝛼<1𝑥,0=0,𝑥∈0,20,𝑡=𝑢2,𝑡=0 (14) with 𝜇=1, 𝐾=1, and 𝑓()=2𝑥,𝑡Γ()𝑥()𝑡2.32−𝑥^1.3+𝑥(2−𝑥)𝑡²+2𝑡². (15)

The exact solution of (14) is 𝑢(𝑥,𝑡)=𝑥(2−𝑥)𝑡². (16)

With 𝜖=10⁻⁶, 𝑇=1.0, 𝑥_𝑅=2.0, and 𝑃=3, the comparison between exact solution and the presented DD algorithm is shown in Table 1. We can find that the DD algorithm compares well with the exact solution.

Table 1 

Comparing exact solution and DD algorithm.

We can replace the DDM (lines 16–27 of Algorithm 1) with Jacobi method. The Jacobi method for a time step has the same parallelism with the DD algorithm. But the Jacobi method needs more iterations. With 𝜖=10⁻⁶ and 𝑃=3, the comparison between Jacobi method and the presented DD algorithm is shown in Table 2. The sum of “count” (total iterations) for all time steps is recorded. We can see that the DDM needs much less iterations than Jacobi method.

Table 2 

Comparing Jacobi method and DDM.

As a part of the future work, we would like to implement an efficient DDM for time fractional equations on parallel computer systems, for example, Tianhe-1A supercomputer [26].

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

This research work is supported by the National Natural Science Foundation of China under Grant no. 11175253. The authors would like to thank the anonymous reviewers for their helpful comments as well.