NUMERICAL DISCRETIZATION AND FAST APPROXIMATION OF A VARIABLY DISTRIBUTED-ORDER FRACTIONAL WAVE EQUATION

We investigate a variably distributed-order time-fractional wave partial differential equation, which could accurately model, e.g., the viscoelastic behavior in vibrations in complex surroundings with uncertainties or strong heterogeneity in the data. A standard composite rectangle formula of mesh size σ is firstly used to discretize the variably distributed-order integral and then the L-1 formula of degree of freedom N is applied for the resulting fractional derivatives. Optimal error estimates of the corresponding fully-discrete finite element method are proved based only on the smoothness assumptions of the data. To maintain the accuracy, setting σ = O(N−1) leads to O(N) operations of evaluating the temporal discretization coefficients. To improve the computational efficiency, we develop a novel time-stepping scheme by expanding the fractional kernel at a fixed fractional order to decouple the fractional operator from the variably distributed-order integral. Only O(log N) terms are needed for the expansion without loss of accuracy, which consequently reduce the computational cost of generating coefficients from O(N) to O(N log N). Optimal-order error estimates of this time-stepping scheme are rigorously proved via novel and different techniques from the standard analysis procedure of the L-1 methods. Numerical experiments are presented to substantiate the theoretical results. Mathematics Subject Classification. 35R11, 65N30. Received March 3, 2021. Accepted August 17, 2021.


Introduction
We begin with a viscoelastic vibratiion model. Consider a homogeneous membrane tautly stretched in a bounded domain Ω in the horizontal plane in its equilibrium position, which undergoes small transverse vibrations in the vertical direction due to the influence of the loading ( , ) or the nonzero initial position or initial velocity. Let ( , ) denote the (positive upward) displacement of the membrane at point ∈ Ω at time in the vertical direction. Because a membrane provides no resistance to bending, the restoration force results only from the tension of the membrane. Let be the constant tension per unit length of the membrane and be the mass per unit area of the membrane. Let ∆ be a small area in the membrane and its projection in the horizontal plane is ∆Ω. Apply Newton's second law to the dynamic equilibrium of vertical forces in the beam Keywords and phrases. Variably distributed-order time-fractional wave equation, viscoelastic problem, well-posedness and regularity, finite element method, optimal-order error estimate, fast algorithm. where refers to the outward normal differential operator to the boundary Ω. Applying Green's formula to the first term on the right-hand side, diving the equation by |∆Ω| and then taking the limit as the diameter of ∆Ω tends to zero gives rise to the classical vibration equation Here := / , := / , and we assume that the membrane is clamped on the boundary Ω. We have assumed that the membrane is clamped on the boundary Ω. However, the classical vibration equation (1.2) does not always provide a physically correct prediction. For instance, model (1.2) predicts that a free vibration (i.e., ≡ 0) will continue forever, which contradicts to the reality that the vibration will die off eventually. The discrepancy is due to the fact that model does not account for the dissipation mechanism in the vibrations. When the membrane is immersed in a viscous medium such as water or air, an additional viscous damping term of the form − ∫︀ ΔΩ , with being the viscous damping coefficient, needs to be included to the right-hand side of equation (1.1) to account for the viscous friction of the medium. On the other hand, if the membrane is attached to elastic material, then the friction assumes a form of − ∫︀ ΔΩ with being the elastic modulus of the medium. Finally, when the membrane is immersed in a viscoelastic medium, such as body fluids in biological and medical fields, then the viscoelastic damping exhibits both viscous property of liquids and elastic character of solids, which shows power-law behavior. Hence, the viscoelastic damping impact can be modeled by − ∫︀ ΔΩ with 0 ≤ ≤ 1, which is well known to accurately describe power-law behaviors [3,27,31]. Here the fractional differential operator := 0 1− for 0 ≤ < 1 with the fractional integral operator 0 1− := ( − /Γ(1 − )) * and := for = 1 [31]. Hence, the viscoelastic damping term − ∫︀ ΔΩ accurately describes the behavior of viscoelastic damping, and includes the elastic resistance and viscoelastic damping as special cases [3,4,8,26,27,31,37,38], and has attracted growing research activities [9,11,13,18,19,22,25,29,33,34,41,42,44]. Consequently, the modeling equation becomes with := / being the viscoelastic damping coefficient and 0 ≤ < 1. When the surrounding viscoelastic medium is highly heterogeneous, the vibration model (1.3) with a viscoelastic damping term of a single fractional order 0 ≤ < 1, which is related to the fractal dimension of the medium via the Hurst index [30], often does not suffice to characterize the impact of the medium. Distributedorder fractional differential operators were introduced to accommodate the integrated effect of the fractional differential operators with respect to a spectrum of fractional orders [6,7,24], which have attracted growing research activities [10,15]. In applications the cyclic motions of the material may change the complex structure of media, leading to a variably distributed-order fractional wave equation [32,36] (1.4) Here Ω ⊂ R ( = 1, 2, 3) is a -dimensional convex polytope, ( ) ≥ 0, ∇ 2 represents the -dimensional Dirichlet Laplacian, and Over the past few decades, there are fruitful progresses on the mathematical and numerical analysis of distributed-order fractional partial differential equations (PDEs) as well as the corresponding fast solvers. For distributed-order problems with = ( ), traditional analysis techniques like the Laplace transform could be used to analyze their well-posedness and regularity [5,7,14,23], which may not be applicable if is time dependent as in the current context. In a recent work the spectral decomposition approach was adopted to analyze a variable-order fractional PDE [43], which circumvents the difficulties caused by the time-dependent variable order and the ideas may work for the propsoed model. Numerical methods and associated analysis for distributed-order PDEs have been extensively investigated [10,21,28]. A typical discretization of distributedorder fractional derivative is to apply the quadrature rules for distributed-order integrals to obtain multi-term fractional derivatives, and then use well-developed schemes like the L-1 methods to discretize each singleorder fractional derivative. However, such discretizations lead to significant computational costs for the case of variably distributed-order fractional derivatives as mentioned in the abstract. Though there are well-developed fast algorithms for distributed-order space-fractional derivatives, see e.g., [17,35], the corresponding studies for variably distributed-order time-fractional problems are meager.
In this paper we develop and analyze a numerical discretization of problem (1.4) and its fast evaluation. In Section 2 we analyze problem (1.4). In Section 3 we develop and analyze a numerical discretization of problem (1.4). In Section 4 we develop and analyze a fast solver of the discretization. In Section 5 we carry out numerical experiments to substantiate the numerical discretization and its fast solver. We finally present concluding remarks for possible extensions of the current work.

Let
(Ω) with ∈ N be the space of continuously differentiable functions of order onΩ, (Ω) with ∈ R ≥0 be the fractional Sobolev space of order on Ω, 0 (Ω) be the subspace of (Ω) enforced with the proper zero boundary conditions, and (ℐ; ), with ℐ being a time interval and being a Banach space, be the space-time space, all equipped with the conventional norms [1,12].

A corresponding fractional wave ordinary differential equation
We prove the well-posedness and regularity analysis of a variably distributed-order fractional wave ordinary differential equation Here > 0 is a prescribed constant and 0 ,ˇ0 and ( ) are given data.
If ∈ 1 [0, ], then ′′′ ∈ (0, ] and for any ∈ (0, ], where the constant is independent of , 0 ,ˇ0, or . Proof. We move ( ) (2.7) From (2.5), we express ′ ( ) in the form and thus the convolution of the third term with cos( ) could be expressed as By setting = ′′ and incorporate the preceding two equations into (2.7) to rewrite (2.7) as a Volterra integral equation as follows (2.8) By the assumptions of the theorem, the kernels of all the right-hand side terms are continuous. We apply the theory of the second kind Volterra integral equation ( [16], Thm. 2.1.1) and use the fact that Consequently, ( ) = ( ) * +ˇ0 + 0 ∈ 2 [0, ] is the solution to (2.2) and all the derivations among (2.6)-(2.8), such as interchanging the order of integration with differentiation or another integration, are applicable. Hence, estimate (2.3) with = 2 holds by (2.9). The uniqueness of the 2 solution to (2.2) follows from that of the continuous solution to the integral equation (2.8).
To estimate ′ , we differentiate (2.6) to get Applying Gronwall's inequality yields (2.3) with = 1. If ∈ 1 [0, ], all terms on the right-hand side of (2.8) are continuously differentiable. This implies ∈ 1 (0, ]. We differentiate (2.8) to obtain By the Condition and the fact Thus the proof is completed.

Finite element approximation and error estimates
To develop and analyze a finite element approximation to problem (1.4), we adopt the order reduction approach [20] to reformulate problem (1.4) in terms of the following first-order system (3.1)
Proof. We use Theorem 2.2 to bound andˆin (3.2) by We may bound ∇ˆsimilarly. We combine the derivative of ( , ) 0 1− ( , ) with respect to and Theorem 2.2 to find that which leads to We then bound by applying Theorem 2.2 Finally, we estimate andˆ1 We could bound and Δˆs imilarly and thus complete the proof.
We follow the ideas in [45] to prove error estimate in the following theorem. , there exists a ∆ 0 > 0 such that the following optimal-order error estimate holds for where is given in Theorem 3.1 and is independent of 0 ,ˇ0, , ∆ , or ℎ.

Remark 3.3.
From the proof we observe that ∆ 0 → 0 if → 1 − , which leads to → ∞. In practice, the fractional derivative of order 0 < < 1 describes the viscoelastic damping and thus we set < 1, which avoids the degeneration of ∆ 0 . The case of = 1 is much more intricate and will be studied in the future work.

A novel discretization scheme and its analysis
In the previous section, we first discretize the variably distributed-order integral in (3.1) by a composite rectangle formula and then discretize each resulting constant-order fractional integral to approximate the variably distributed-order operator, and the corresponding finite element scheme (3.8)-(3.9) has an optimal-order error estimate ( + ∆ + ℎ 2 ). Numerically, we need to compute all the coefficients , for 1 ≤ ≤ ≤ and 1 ≤ ≤ . To maintain the first-order accuracy in time, we may set = (∆ ), which leads to ( 3 ) computations to generate the coefficients { , } that is computationally consuming. Motivated by this observation, we consider a novel time-stepping scheme by discretizing the variably distributed-order fractional integral via a different approach. We first discretize the fractional-order integral at = by where the coefficients˜, and local truncation error˜′ are given bỹ .

(4.2)
The key is to develop a fast approximation method of {˜, }. Define a function where the local truncation error * , is given by Therefore,˜, with − ≥ 2 in (4.2) can be approximated by˜, defined bỹ and the corresponding error ** , :=˜, −˜, ( − ≥ 2) is given by Substituting˜, by˜, in (4.1) we obtain a novel approximation of variably distributed-order fractional integral with the local truncation error˜: We multiply the governing equation by ∈ 1 0 (Ω) on Ω to obtain that for any ∈ 1 0 and = 1, 2, · · · , holds. Then the following estimates hold where is given in Theorem 3.1 and is a constant independent from data.
Combining thsi with the definition of ** , in (4.4) and the Stirling's formula ( + 1)! > ( + 1) +3/2 −( +1) we find (4.11) Thus in order to retain the first-order accuracy in time, we require which is equivalent to Let + 1 = ln for some > 1 and thus Therefore it suffices to select satisfying ( − 1) ≥ 2, and (4.11) reduces to ‖˜′ ′ ‖^∞ ( 2 ) ≤ ∆ and consequently Similar to (4.5), we decompose˜1 (4.12) The first two terms on the right-hand side can be estimated by By Condition and similar techniques in (4.11), we bound the third term on the right-hand side of (4.12) by Therefore, we finish the proof of the second statement of (4.10). where is a positive constant independent from data.
Proof. By (4.5) with ≡ 1 we have By Condition S and a similar derivation as (4.11) we have (4.14) Similar to the proof of (3.16), we obtain ∑︀ * =˜, ≤ . The second term on the right-hand side of (4.14) can be bounded via Thus we obtain from (4.14) that ≤ , which completes the proof. . Then there exists a ∆ 0 > 0 such that the optimal-order error estimate holds for 0 < ∆ ≤ ∆ 0

Error estimates and computational efficiency of the novel time-stepping scheme
where is given in Theorem 3.1 and is independent of 0 ,ˇ0, , ∆ or ℎ.
We summarize the above observations to complete the proof.
Remark 4.5. Based on the proof of Theorem 4.4, it suffices to approximate˜, by˜, for − ≥ rather than − ≥ 2, with the same magnitude of computational cost.

Numerical experiments
We carry out numerical experiments to investigate the performance of the finite element method (FEM) In numerical experiments, we apply the uniform rectangular partition on Ω with the mesh size ℎ, and set = ∆ and = ⌊ 3/2 ln ⌋ = −⌊ 3/2 ln ∆ ⌋ for the FEM and nFEM, respectively. As the spatial discretization is standard, we set ℎ = 2 −6 and only measure the temporal convergence rates by and similarly, ‖ −˜‖^∞ (^2) and ‖ −˜‖^∞ (^2) with respect to˜and˜, respectively. We also record the following two kinds of CPU times of generating entries of temporal discretization coefficients, that is, let and denote the CPU time of computing coefficients ofˆ1 for 1 ≤ ≤ , respectively. The increasing ratio of CPU times is defined by where is the number of degree of freedom of the temporal partition in the th computation and is the corresponding order of magnitude.

Model (1.4) with smooth solutions
We test the performance of FEM and nFEM with smooth solutions given by ( ). ( , ) = − sin(2 ) sin( ), ( ). ( , , ) = − sin(2 ) sin( ) sin( ), and the right-hand side terms ( , ) or ( , , ) are computed accordingly. We present errors of numerical solutions in Tables 1-2, from which we observe that both FEM and nFEM have the same accuracy. We also present and for case (i) in Table 3, from which we observe that nFEM is more efficient that FEM. We further plot these CPU times in Figure 1, which indicates that the increases cubically while the increases almost quadratically, which is highly consistent with the analysis in Theorem 4.4.
As the exact solutions are not available, we set the numerical solutions of FEM with = 2 10 and ℎ = 2 −6 to be the reference solutions. The errors and convergence rates of the schemes are presented in Tables 4-5, which again demonstrate the accuracy of the nFEM.   Table 3. CPU times of generating temporal discretization coefficients.

Concluding remarks
In this paper we prove the well-posedness and regularity of a variably distributed-order time-fractional vibration PDE (1.4), and develop and analyze a fast numerical discretization of the model. In the mathematical model the fractional order time derivative term accurately describes the power-law behavior of viscoelastic damping mechanisms in the vibrations. The time-dependent variably distributed-order time-fractional derivative term accounts for the integrated effect of the fractional differential operators with respect to a spectrum of fractional orders and its impact on the evoluation of the viscoelastic properties of the materials.
In the vibration model (1.4) the viscoelastic damping coefficient and the variably distributed order are assumed to be spatially homogeneous based on the following motivations: The model (1.4) describes viscoelastic vibrations of physical structures such as membranes (and strings), which have negligible thickness (and width) and so provide only negligible resistance to bendings. In these circumstances the variation in thickness (and width) in membrances (and strings) shall not introduce spatial variations of and . Hence, the spatial variations of the these parameters probably just introduce extra complexities in terms of measurements of the  corresponding parameters and the design and manufacture of the physical structures, and the model (1.4) suffices in many applications. Mathematically and numerically, these simplifications simplifies the analysis of the model and its numerical discretizations. Space-dependent parameters do occur in complicated scenarios. For instance, long time vibrations may lead to material fatigue or even damage such as locally reduced stiffness or tension, which may yield spatial heterogeneities. The corresponding mathematical model may assume the form of (1.4) where and are also space dependent, with possible nonlinearities and couplings to other equation or constraints. Additional issues that need to be addressed include the following: (i) The spectral decomposition based mathematical analysis in this paper does not apply to this case. (ii) Numerically, due to the coupling of , and the inner product in the FEM, the corresponding temporal discretization coefficients could not be separated from the inner product, which could complicate the numerical analysis. Finally, it is possible to combine the novel discretization and the fast algorithm developed and analyzed in this paper, which are devoted only for the evaluation of the variably distributed-order time-fractional derivative, with the existing numerical techniques on complex domains to analyze the proposed model on a more complex geometry. All the possible extensions will be investigated in the near future.