• Views 745
• Citations 3
• ePub 16
• PDF 506
Mathematical Problems in Engineering
Volume 2014, Article ID 194249, 5 pages
http://dx.doi.org/10.1155/2014/194249
Research Article

## An Efficient Iteration Method for Toeplitz-Plus-Band Triangular Systems Generated from Fractional Ordinary Differential Equation

1College of Aerospace Science and Engineering, National University of Defense Technology, Changsha 410073, China
2Science and Technology on Space Physics Laboratory, Beijing 100076, China
3School of Computer Science, National University of Defense Technology, Changsha 410073, China

Received 2 January 2014; Accepted 27 May 2014; Published 12 June 2014

Copyright © 2014 Chunye Gong et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### Abstract

It is time consuming to numerically solve fractional differential equations. The fractional ordinary differential equations may produce Toeplitz-plus-band triangular systems. An efficient iteration method for Toeplitz-plus-band triangular systems is presented with computational complexity and memory complexity in this paper, compared with the regular solution with computational complexity and memory complexity. is the discrete grid points. Some methods such as matrix splitting, FFT, compress memory storage and adjustable matrix bandwidth are used in the presented solution. The experimental results show that the presented method compares well with the exact solution and is 4.25 times faster than the regular solution.

#### 1. Introduction

Fractional differential equation (FDE) plays an important role in dynamical systems [1] and has more than 300 years of research history [2]. Many analytical solutions and numerical solutions [36] have been proposed for FDE, such as finite difference method [7, 8], finite element method [9], and spectral method [10, 11]. In recent times, interest in fractional ordinary differential equations (FODE) has increased [1215]. The derivatives in the FODE are approximated by linear combinations of function values at the discrete grid points. Compared with integer ordinary differential equations, the FODE has nonlocal effect, which means a grid point may rely on the grid points far away from its position. And a grid point of the classical integer equations may only rely on its several neighboring grid points.

For integer order equations, the coefficient matrices are often sparse. Because of the nonlocal property of fractional differential operators, the numerical methods for fractional diffusion equations often generate dense or even full coefficient matrices [16]. This nonlocal property makes the computation of FODE and FDE much heavier than that of the traditional integer equations. The short memory principle [17], parallel computing [1821], fast Fourier transformation (FFT) [22, 23], multigrid method [24], and preconditioner technologies [25, 26] are used to overcome this heavy computation. Gong et al. presented many parallel algorithms for different FDEs on both traditional and heterogeneous parallel platforms [16, 18]. Diethelm [19] proposed a parallel second-order Adams-Bashforth-Moulton method for a FODE. Wang and Du [26] proposed a superfast-preconditioned iterative method for steady-state two-side space-fractional diffusion equations.

The fractional ordinary differential equations may produce Toeplitz-plus-band triangular systems. Toeplitz-plus-band systems were studied by professors Chan and Ng  [27]. They considered the solutions of Hermitian Toeplitz-plus-band systems , where are -by- Toeplitz matrices and are -by- band matrices with bandwidth independent of . and are both Hermitian matrix. The authors proved that if is generated by a nonnegative piecewise continuous function and is positive semidefinite, then there exists a band matrix , with bandwidth independent of , such that the spectra of are uniformly bounded by a constant independent of . The band preconditioner was developed for Hermitian Toeplitz systems [28]. The recursive blocked algorithms were proposed for triangular systems and the recursive algorithms lead to an automatic variable blocking that has the potential of matching the memory hierarchies of today’s HPC (high performance computing) systems  [29, 30].

This paper focuses on the fractional ordinary differential equation [13]: where , , , and . The fractional derivative is in the Caputo form [31].

Define for , where is a positive integer, and are step size. Assume to be the numerical approximation to and the numerical approximation to . Using the Grünwald approximation, the finite difference scheme for (1) is shown as follows: where the normalized Grünwald weight is defined by

Equation (2) results in a linear system of equations where and . If , the term should be included in . is the coefficient matrix. is defined by

#### 2. Method

##### 2.1. Analysis

In a more explicit format, matrix can be represented as

The linear system (4) can be solved with computational complexity , shown in Algorithm 1. The output equals .

Algorithm 1: Forward substitution for lower triangular matrix.

From (6), we can see that has some properties.(1)is a low triangular, diagonal dominant matrix.(2)One has for , . This property is determined by the normalized Grünwald weight and is the mathematical background of short memory principle. This property means that for grid point , if the distance of grid point is smaller than that of grid point , has more impact on than .(3)If is split into two matrices and , . is a banded matrix and the bandwidth (number of diagonals) . Matrix can be factorized into a product . is a diagonal matrix . is a Toeplitz matrix with nonzero diagonals on its left-bottom part.(4)The Toeplitz matrix can be stored with memory space compared with for a general low triangular matrix with order .

##### 2.2. Efficient Iteration Method

Equation (4) evolves as follows:

So the linear algebra can be solved iteratively, shown in (10). Because is a diagonal matrix, keeps associative law and commutative law for matrix-matrix multiplication. The rate of convergence associated with (10) depends on the eigenvalues of the iteration matrix [32]:

Assuming error with satisfies , then with norm [32]:

So the spectral radius of () determines the asymptotic behavior of . From Theorem of [32], we can conclude that if and only if , (10) will converge to . Generally speaking, the iteration is expected to work well with small .

Assume the bandwidth of matrix is and . Solving needs about arithmetical operations. If is near , there are about arithmetical operations with forward substitution. So the computational complexity of is .

Assume , , and . The computation of and needs multiplications and additions, respectively. Because only has nonzero diagonals on its left-bottom part, only the front elements of are effective for the multiplication . The back elements of are zero. So can be regarded as a Toeplitz matrix vector multiplication with order . It is well known that Toeplitz matrix vector multiplication with order can be finished with operations [33]. The Toeplitz matrix vector multiplication can be computed by FFTs by first embedding into a -by- circulant matrix. The cost of circulant matrix vector multiplication is by using FFTs of length .

So the cost of each iteration of (10) is . If is a diagonal dominant matrix, we can expect (10) can converge with not too many iterations. The efficient iteration method is shown in Algorithm 2.

Algorithm 2: The efficient iteration method for FODE.

In Algorithm 2, stands for the value of previous iteration and stands for the current iteration. stands for with . equals the reciprocal of with . stands for with . and are arrays. equals .

The value of can affect the performance of Algorithm 2 shown in Table 1.

Table 1: Impact of .

Algorithm 2 has five features/advantages compared to Algorithm 1.(1)Split the coefficient matrix and solve the triangular system iteratively.(2)Use FFT to compute matrix vector multiplication.(3)Precompute .(4)Compress storage.(5)Adjust parameter .

#### 3. Numerical Example

The experiment platform is a laptop with Intel(R) Core (TM) i3-3110M CPU, 2 GB main memory, and Windows 7 operating system. The CPU clock frequency is 2.40 GHz. The code is developed with MATLAB R2012a and runs on default double precision floating point operations.

The following fractional () ordinary differential equation [13] was considered: where .

The exact solution of (14) is

The efficient iteration method of Algorithm 2 compares well with the exact solution to the FODE in the test case of (14), shown in Figure 1. The is . The maximum absolute error is . The difference between the efficient iteration method and the forward substitution Algorithm 1 is only . The efficient iteration method and the regular forward substitution solution have no noticeable artifacts.

Figure 1: Comparison of exact solution to the solution of the fast solution at time = 1.0.

The performance comparison between regular forward substitution solution of Algorithm 1 and efficient iteration method of Algorithm 2 is shown in Table 2. Columns 2 and 3 of Table 2 are the runtime and the runtime is recorded in seconds. With , the maximum speedup is 4.25. Because the speedup increases with , the bigger is, the higher the speedup that can be expected is. Because of the 2 GB memory limitation, the compress memory usage is also used in Algorithm 1.

Table 2: Performance comparison between regular solution and the presented efficient iteration method.

The impact of on the performance of Algorithm 2 is shown in Table 3. The runtime of the presented method varies with . So is a key parameter for the performance of Algorithm 2. In real fractional ordinary applications, the proper should be chosen.

Table 3: Impact of for .

The presented iteration method should be regarded as an iteration method to solve not only the system generated from FODE but also the more general Toeplitz-plus-band triangular systems. The technology of parallel computing is very useful, but with less mathematical background. Parallel computing is attractive for fractional differential equations [34]. As a part of future work, first, we would like to parallelize the presented solution on shared memory or distributed memory systems. Second, accelerating the presented efficient iteration method on heterogeneous architecture [3538] should also be interesting.

#### 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 in part by the National Natural Science Foundation of China under Grant nos. 11175253 and 61170083, in part by the Specialized Research Fund for the Doctoral Program of Higher Education under Grant no. 20114307110001, and in part by the 973 Program of China under Grant no. 61312701001. The authors would like to thank the anonymous reviewers for their helpful comments as well.

#### References

1. R. Caponetto, G. Dongola, F. Pappalardo, and V. Tomasello, “Auto-tuning and fractional order controller implementation on hardware in the loop system,” Journal of Optimization Theory and Applications, vol. 156, no. 1, pp. 141–152, 2013.
2. J. A. Tenreiro Machado, V. Kiryakova, and F. Mainardi, “A poster about the recent history of fractional calculus,” Fractional Calculus & Applied Analysis, vol. 13, no. 3, pp. 329–334, 2010.
3. A. Atangana and A. Kilicman, “Analytical solutions of the space-time fractional derivative of advection dispersion equation,” Mathematical Problems in Engineering, vol. 2013, Article ID 853127, 9 pages, 2013.
4. S. Yu. Lukashchuk, “An approximate solution method for ordinary fractional differential equations with the Riemann-Liouville fractional derivatives,” Communications in Nonlinear Science and Numerical Simulation, vol. 19, no. 2, pp. 390–400, 2014.
5. D. Baleanu, K. Diethelm, E. Scalas, and J. J. Trujillo, Fractional Calculus: Models and Numerical Methods, vol. 3, World Scientific Publishing, 2012.
6. E. Gao, S. Song, and X. Zhang, “A numerical algorithm for solving a four-point nonlinear fractional integro-differential equations,” Mathematical Problems in Engineering, vol. 2012, Article ID 329575, 11 pages, 2012.
7. S. Zhai, X. Feng, and Z. Weng, “New high-order compact ADI algorithms for 3D nonlinear time-fractional convection-diffusion equation,” Mathematical Problems in Engineering, vol. 2013, Article ID 246025, 11 pages, 2013.
8. J. Huang, Y. Tang, L. Vázquez, and J. Yang, “Two finite difference schemes for time fractional diffusion-wave equation,” Numerical Algorithms, vol. 64, no. 4, pp. 707–720, 2013.
9. N. J. Ford, J. Xiao, and Y. Yan, “A finite element method for time fractional partial differential equations,” Fractional Calculus and Applied Analysis, vol. 14, no. 3, pp. 454–474, 2011.
10. M. Zayernouri and G. E. Karniadakis, “Exponentially accurate spectral and spectral element methods for fractional ODEs,” Journal of Computational Physics, vol. 257, pp. 460–480, 2014.
11. J. Huang, N. Nie, and Y. Tang, “A second order finite difference-spectral method for space fractional diffusion equations,” Science China: Mathematics, vol. 57, no. 6, pp. 1303–1317, 2014.
12. J. Leszczyński and M. Ciesielski, “A numerical method for solution of ordinary differential equations of fractional order,” in Parallel Processing and Applied Mathematics, pp. 695–702, Springer, 2006.
13. X. Cai and J. H. Chen, “Numerical method for fractional order ordinary differential equation,” Journal of Jimei University(Natural Science), vol. 12, no. 4, pp. 367–370, 2007.
14. R. Lin and F. Liu, “Fractional high order methods for the nonlinear fractional ordinary differential equation,” Nonlinear Analysis: Theory, Methods & Applications, vol. 66, no. 4, pp. 856–869, 2007.
15. J. Cao and C. Xu, “A high order schema for the numerical solution of the fractional ordinary differential equations,” Journal of Computational Physics, vol. 238, pp. 154–168, 2013.
16. C. Gong, W. Bao, and G. Tang, “A parallel algorithm for the Riesz fractional reaction-diffusion equation with explicit finite difference method,” Fractional Calculus and Applied Analysis, vol. 16, no. 3, pp. 654–669, 2013.
17. W. Deng, “Short memory principle and a predictor-corrector approach for fractional differential equations,” Journal of Computational and Applied Mathematics, vol. 206, no. 1, pp. 174–188, 2007.
18. J. Liu, C. Gong, W. Bao, G. Tang, and Y. Jiang, “Solving the caputo fractional reaction-diffusion equation on GPU,” Discrete Dynamics in Nature and Society. In press.
19. K. Diethelm, “An efficient parallel algorithm for the numerical solution of fractional differential equations,” Fractional Calculus and Applied Analysis, vol. 14, no. 3, pp. 475–490, 2011.
20. C. Gong, W. Bao, G. Tang, Y. Jiang, and J. Liu, “A parallel algorithm for the two dimensional time fractional diffusion equation with implicit difference method,” The Scientific World Journal, vol. 2014, Article ID 219580, 8 pages, 2014.
21. C. Gong, W. Bao, G. Tang, Y. Jiang, and J. Liu, “A domain decomposition method for time fractional reaction-diffusion equation,” The Scientific World Journal, vol. 2014, Article ID 681707, 5 pages, 2014.
22. H. Wang and K. Wang, “An $O\left(N{log}^{2}N\right)$ alternating-direction finite difference method for two-dimensional fractional diffusion equations,” Journal of Computational Physics, vol. 230, no. 21, pp. 7830–7839, 2011.
23. H. Wang, K. Wang, and T. Sircar, “A direct $O\left(N{log}^{2}N\right)$ finite difference method for fractional diffusion equations,” Journal of Computational Physics, vol. 229, no. 21, pp. 8095–8104, 2010.
24. H.-K. Pang and H.-W. Sun, “Multigrid method for fractional diffusion equations,” Journal of Computational Physics, vol. 231, no. 2, pp. 693–703, 2012.
25. S.-L. Lei and H.-W. Sun, “A circulant preconditioner for fractional diffusion equations,” Journal of Computational Physics, vol. 242, pp. 715–725, 2013.
26. H. Wang and N. Du, “A superfast-preconditioned iterative method for steady-state space-fractional diffusion equations,” Journal of Computational Physics, vol. 240, pp. 49–57, 2013.
27. R. H. Chan and K.-P. Ng, “Fast iterative solvers for Toeplitz-plus-band systems,” SIAM Journal on Scientific Computing, vol. 14, no. 5, pp. 1013–1019, 1993.
28. R. H. Chan and P. T. P. Tang, “Fast band-Toeplitz preconditioners for Hermitian Toeplitz systems,” SIAM Journal on Scientific Computing, vol. 15, no. 1, pp. 164–171, 1994.
29. I. Jonsson and B. Kågström, “Recursive blocked algorithm for solving triangular systems. I. One-sided and coupled Sylvester-type matrix equations,” ACM Transactions on Mathematical Software, vol. 28, no. 4, pp. 392–415, 2002.
30. I. Jonsson and B. Kågström, “Recursive blocked algorithm for solving triangular systems. II. Two-sided and generalized Sylvester and Lyapunov matrix equations,” ACM Transactions on Mathematical Software, vol. 28, no. 4, pp. 416–435, 2002.
31. I. Podlubny, Fractional Differential Equations, vol. 198, Academic Press, San Diego, Calif, USA, 1999.
32. G. H. Golub and C. F. Van Loan, Matrix Computations, vol. 3, Johns Hopkins University Press, Baltimore, Md, USA, 4th edition, 2013.
33. R. H.-F. Chan and X.-Q. Jin, An introduction to iterative Toeplitz solvers, vol. 5, Society for Industrial and Applied Mathematics (SIAM), Philadelphia, Pa, USA, 2007.
34. C. Gong, W. Bao, G. Tang, B. Yang, and J. Liu, “An effcient parallel solution for Caputo fractional reaction-diffusion equation,” The Journal of Supercomputing, 2014.
35. X. J. Yang, X. K. Liao, K. Lu, Q. F. Hu, J. Q. Song, and J. S. Su, “The TianHe-1A supercomputer: its hardware and software,” Journal of Computer Science and Technology, vol. 26, no. 3, pp. 344–351, 2011.
36. Q. Wu, C. Yang, T. Tang, and L. Xiao, “Exploiting hierarchy parallelism for molecular dynamics on a petascale hetero geneous system,” Journal of Parallel and Distributed Computing, vol. 73, no. 12, pp. 1592–1604, 2013.
37. C. Gong, J. Liu, L. Chi, H. Huang, J. Fang, and Z. Gong, “GPU accelerated simulations of 3D deterministic particle transport using discrete ordinates method,” Journal of Computational Physics, vol. 230, no. 15, pp. 6010–6022, 2011.
38. C. Gong, J. Liu, H. Huang, and Z. Gong, “Particle transport with unstructured grid on GPU,” Computer Physics Communications, vol. 183, no. 3, pp. 588–593, 2012.