A fast multiscale Galerkin algorithm for solving boundary value problem of the fractional Bagley–Torvik equation

In this paper, a fast multiscale Galerkin algorithm is developed for solving the boundary value problem of the fractional Bagley–Torvik equation. For this purpose, we employ multiscale orthogonal functions having vanishing moments as the basis of the trial space, and we propose a truncation strategy for the coefficient matrix of the corresponding discrete system which leads to a fast algorithm. We show the algorithm has nearly linear computational complexity (up to a logarithmic factor). Numerical experiments are presented to illustrate the efficiency, accuracy and convergence of the proposed algorithm. Also, comparisons with some other existing methods are made to confirm the reliability of the algorithm

In this paper, we consider the following Bagley-Torvik equation: subject to boundary conditions: where θ , σ and 0 < α < 1 are constants, and f (t) is a continuous function defined on the interval [0, 1]. Here the notation D α represents the Caputo fractional derivative defined by where Γ (·) is the gamma function defined by The Bagley-Torvik equation was first introduced by Bagley and Torvik and treated as a mathematical model of the motion of a thin rigid plate immersed in a Newtonian fluid [22]. It was then widely used to simulate the dynamic responses of viscoelastically damped structures [23]. Many numerical approaches have been developed for solving the Bagley-Torvik equation. These methods include the modified Galerkin method [10], collocation methods [12,13], the wavelet method [14,15], the finite difference method [16], the spline method [17][18][19] and the hybridizable discontinuous Galerkin method [20], and the operational Tau method [21]. However, due to the non-local property of the fractional differential operators, the numerical methods for fractional equation including Bagley-Torvik equation lead to dense coefficient matrices. When the dimension of coefficient matrix is large, the computational cost for generating the matrix and then solving the corresponding linear system is huge.
The main purpose of this paper is to develop a fast multiscale Galerkin algorithm for solving the Bagley-Torvik equation. We introduce a matrix truncation strategy by choosing multiscale orthogonal basis functions having vanishing moments, which forms a basis for the fast algorithm for solving Eq. (1). Specifically, the multiscale orthonormal basis constructed in [24] is employed to discretize the Bagley-Torvik equation, which leads to a linear system with a numerically sparse coefficient matrix. That is, the absolute values of most of the entries of the matrix are very small. We then set the entries with small value to zero by a matrix truncation strategy, which yields a sparse matrix. We show that the number of nonzero entries of the truncated matrix is linear (up to a logarithmic factor) with respect to the dimensions of the matrix. This method cannot only generate the coefficient matrix rapidly, but also make it easy to solve the resulting linear system with sparse coefficient matrix. This paper is organized as follows. In Sect. 2, the multiscale orthonormal bases in Sobolev space H 1 0 (0, 1) is introduced. In Sect. 3, the fast multiscale Galerkin algorithm with a practical matrix truncation strategy is proposed to solve the Bagley-Torvik equation boundary value problem. Some numerical examples are presented in Sect. 4, and the conclusion is included in the last section.

Multiscale orthonormal bases
The multiscale orthonomal bases for Sobolev space on the unit interval I = [0, 1] were constructed by Chen, Wu and Xu in Ref. [24] for solving differential equations. For the convenience of the reader, we briefly introduce the construction and properties of this multiscale orthonomal bases, which are required necessarily for designing the fast algorithm for the Bagley-Torvik equation.
We start with some useful notations. Let N := {1, 2, . . . }, N 0 := N ∪ {0} and Z n := {0, 1, 2, . . . , n -1}, for n ∈ N. We denote by (·, ·) the inner product on the space L 2 (I) with the L 2 norm · 2 . Let H 1 0 (I) denote the Sobolev space of elements u vanishing at the endpoints that u(0) = u(1) = 0. The inner product and norm of H 1 0 (I) are defined by and respectively. Let n be the uniform mesh which divides the interval I into 2 n pieces, and let X n be the piecewise polynomial space of order k associated with the partition n . It is concluded that the sequence of X n is nested, i.e., which yields the decomposition where W n is the orthogonal complement of X n-1 in X n in the sense of the inner product ·, · 1 . Repeating this decomposition, we have the multiscale space decomposition as follows: where W 0 := X 0 . Denoting x(n) := dim X n and w(n) := dim W n , by the definition of X n we have The following lemma shows that when the bases of X 0 and W 1 are constructed, the basis of the space X n can be recursively generated, and the sequence of X n is ultimately dense in the Sobolev space H 1 0 (I).
Then the orthonormal basis of the space W i+1 can be recursively constructed by Furthermore, We define J n := {(i, j), j ∈ Z w(i) , i ∈ Z n+1 } and assume that W i has an orthonormal basis {w ij : j ∈ Z w(i) }, then {w ij : (i, j) ∈ J n } forms an orthonormal basis of X n . The basis has several important properties, we describe parts of them that are relevant to our algorithm.
(P1) (Compact support) The support S ij := supp w ij holds that where δ ii is the Kronecker delta. (P3) (Vanishing moment) For any i ≥ 1 and p ∈ Z k-1 ,

This, with integration by parts, implies that
To end this section, we give the concrete multiscale linear basis functions which will be used in our numerical experiments. For linear case, k = 2, dim(W i ) = 2 i-1 , for i ≥ 1. Obviously, 0 is the only linear function which vanishes at both 0 and 1. That is X 0 = {0}. The basis of the space W 1 is given by (see [24]) Using w 10 (t) and Lemma 2.1, we can generate the bases of the spaces W 2 , W 3 , . . . , W n in turn, and then we obtain the basis of the approximation subspace X n . We illustrate in Fig. 1 the graph of the basis functions for X 6 .

Fast multiscale Galerkin algorithm
In this section, we first present a Galerkin scheme associate with the multiscale basis introduced in last section for solving Eq. (1). Then a matrix truncation strategy is proposed, which leads to the fast multiscale Galerkin algorithm. We begin with the variational (weak) formulation of the (1) subject to the boundary conditions u(0) = u(1) = 0, which reads: find u ∈ H 1 0 (I) such that Making use of the basis of X n , the Galerkin scheme for solving Eq. (2) is to seek a vector To distinguish the method from the traditional Galerkin method, the scheme (3) is called the multiscale Galerkin method due to the use of the multiscale basis {w ij , (i, j) ∈ J n }. Let From the properties (P1) and (P2), it is easily seen that both matrices I n and E n are sparse. In fact, I n is an identity. However, D n is a dense matrix because of the non-local property of the fractional order differential operator. Moreover, all the entries of D n are defined by singular integrals Therefore, it is computationally costly to establish such a matrix numerically, which becomes an obvious computational deficiency of the Galerkin method for solving such a kind of equations. Fortunately, we observe that the matrix D n is numerically sparse under the multiscale orthonormal basis. That is, a large number of the entries are very small in  Fig. 2 the matrix D n with respect to the multiscale basis with n = 7 and α = 1 2 . According to the different scales of spaces, we partition matrix D n into a block matrix We can see from Fig. 2 that the absolute values of the entries lying outside the diagonals of the blocks are small, which is similar to the case of multiscale methods for solving the Fredholm integral equation [25][26][27]. This observation motivates us to "truncate" the small entries to zero by the same matrix truncation strategy as the references [25][26][27]. Specifically, For (i, j), (i , j ) ∈ J n , we defineD where dist(S ij , S i j ) is the distance between S ij and S i j , and the truncation parameter ε n ii are chosen as [25] ε n ii := max μ2 -n+λ(n-i)+λ (n-i ) , ρ 2 -i + 2 -i for some positive constants μ, λ, λ and ρ > 1. Then we obtain a truncation matrix Replacing D n in Eq. (4) byD n leads to the fast multiscale Galerkin algorithm (FMGA): find c n := [c ij , (i, j) ∈ J n ] such that Solvingc n from Eq. (6), we then obtain the fast multiscale Galerkin solutioñ u n (t) := (i,j)∈J nc ij w ij (t).
Similar to the analysis of [26,27] or Theorem 3.3 of [25], we have the following lemma about the computational complexity of the fast multiscale Galerkin scheme (6), which is measured by the number of nonzero entries of the truncation matrixD n .

Theorem 3.1 For any n ∈ N, we have
where N (M) denotes the number of nonzero entries of matrix M.
The above theorem shows that the truncation matrixD n is sparse, which is very beneficial to the solution of the resulting large linear system by iterative method. More importantly, it greatly reduces the computational burden for calculating all the entries of matrix, in other words, it greatly saves the time for generating the coefficient matrix.

Numerical examples
In this section, we present some numerical illustrations for the solution of the boundary value problem of the Bagley-Torvik equation to show the accuracy and the efficiency of the proposed method. All the numerical calculations are implemented with Matlab 2012a in OS Windows 10 with 2.5 G CPU and 8 G memory.
Example 1 Consider the following Bagley-Torvik equation [17][18][19]: The exact solution of this equation is u * (t) = t λ-1 (t -1). We use the linear multiscale basis to discretize the equation. In order to illustrate the computational efficiency of the FMGA, we solve the equation by both the original scheme (4) (OMGM) and the truncated scheme (6) (FMGA). In all three examples, we choose μ = ρ = 2, λ = 1 and λ = 5 6 in the truncation parameter ε n ii . The numerical results for θ = 0.5, σ = 1, λ = 5 and α = 0.5 are reported in Tables 1 and 2. In the two tables, n, x(n) and T (T) denote the level number of approximation space, the dimension of approximation space and the time of OMGM (FMGA) for generating the coefficient matrix, respectively. The notations · 1 and · 2 stand for the H 1 0 -error and L 2 -error, respectively, and the corresponding convergence orders (C.O.) are computed by the formulas C.O. := log 2 u * -ũ n-1 p u * -ũ n p , C.O.:=log 2 u *u n-1 p u *u n p , where p = 1 or 2;ũ n and u n are the numerical solutions solved by FMGA and OMGM, respectively. It is seen from Table 1 that FMGA and OMGM have the same optimal convergence orders 1 and 2 in the H 1 0 norm and L 2 norm, respectively, and they have almost the same accuracy, which implies that our truncation strategy does not affect the overall accuracy. However, the computational time in Table 2 and Fig. 3 indicate that FMGA is remarkably faster than OMGM. To measure the degree of sparsity of the truncation matrix, Table 1 Numerical results for Example 1 (θ = 0.5, σ = 1, λ = 5, α = 0.5)   we denote by C.R. in Table 2 the ratio of number of nonzero entries to the total number of entries ofD n . We observe that the larger the dimension ofD n , the smaller the ratio. For example, C.R. is only 0.0389 or 3.89% for n = 10, which means you can only calculate less than 4 out of 100 entries. To visualize the sparseness of the matrix, we plot in Fig. 4 the distribution of nonzero entries of D 10 andD 10 , in which the nonzero entries are identified by a black region. A comparison of the absolute error between our method (FMGA) and the methods in [17][18][19] is reported in Table 3, which reveals that the our results have superiority in accuracy.
Example 2 Consider the following Bagley-Torvik equation [17,28]:  Table 3 Comparison of the absolute error with methods in [17][18][19] for Example 1 x Our method GD in [17] NSGDin [17] 3-WSGDin [17] Methodin [18] Methodin [19] 0  Example 3 As the last example we consider the equation [12,17,28] The exact solution is The numerical results obtained by FMGA are presented to show the accuracy, convergence and computing time in Table 4, Fig. 6 and Fig. 7. We also compare the absolute error of our method with that of the methods in [12,17,28] and listed in Table 5, where n * -1 and m are the number of the spline and Bessel functions used in the approximate solution in [17] and [12,28], respectively. These results show that our algorithm has better accuracy with fewer basis functions (x(3) = 7).

Conclusion
In this paper, a fast multiscale Galerkin algorithm based on a matrix truncation strategy for solving boundary value problem of the Caputo fractional Bagley-Torvik equation is presented. The algorithm reduces the computational complexity from O(N 2 ) to O(N log N) (N is the dimension of the approximation space), which makes the algorithm efficient. Examples are presented to illustrate the performance of our algorithm by comparing it with the original Galerkin scheme (OMGM) from accuracy, convergence order and computing time. It shows that they have the same convergence orders and accuracy, but FMGA is much faster than OMGM. Numerical results are also compared with some recent methods, which demonstrate the effectiveness of our algorithm and its superiority in accuracy. The proposed algorithm still needs rigorous theoretical analysis, which constitutes the subject of our ongoing work.