Fourth-order spatio-temporal difference scheme for the Riesz space distributed-order advection–dispersion equation

1 Introduction

This article presents a numerical investigation of the Riesz space distributed-order advection–dispersion equation [1].where p ( α ) and q ( β ) are, respectively, non-negative and bounded weight functions that satisfy the conditions

∂ u x , t ∂ t = ∫ 0 1 ⁢ p ⁡ α ∂ α u x , t ∂ | x | α ⁢ d α + ∫ 1 2 ⁢ q ⁡ β ∂ β u x , t ∂ | x | β ⁢ d β + f ⁡ x , t , x , t ∈ a , b × 0 , T , u ⁡ x , 0 = φ ⁡ x , x ∈ a , b , u ⁡ a , t = 0 , ⁢ u ⁡ b , t = 0 , t ∈ 0 , T , (1)

p α ≥ 0 , p α ≢ 0 , α ∈ 0,1 , 0 < ∫ 0 1 p α d α < ∞ , q β ≥ 0 , q β ≢ 0 , β ∈ 1,2 , 0 < ∫ 1 2 q β d β < ∞ . .

The Riesz fractional derivative∂ ζ u ( x , t ) ∂ | x | ζ ( ζ = α , β ) is given by [2]:where D x , − ζ u ( x , t ) and D x , + ζ u ( x , t ) are the Riemann–Liouville (RL) fractional derivatives [3]. These are defined as

∂ ζ u x , t ∂ | x | ζ = − 1 2 ⁡ cos ζ π 2 D x , − ζ u x , t + D x , + ζ u x , t ,

D x , − ζ u x , t = 1 Γ n − ζ ∂ n ∂ x n ∫ a x u τ , t x − τ ζ − 1 d τ , D x , + ζ u x , t = − 1 n Γ n − ζ ∂ n ∂ x n ∫ x b u τ , t τ − x ζ − 1 d τ , n = ⌈ ζ ⌉ .

Over the past few decades, considerable attention has been paid to anomalous transport, solute migration in heterogeneous media, and viscoelastic dissipation observed across diverse fields, including the natural sciences, biology, geosciences, medicine, signal processing, and finance [3, 4, 5, 6, 7]. Because classical Gaussian statistics are violated in these processes, traditional integer-order advection–dispersion equations fail to provide accurate descriptions [1, 2]. To overcome this drawback, single-term and multi-term fractional differential equations, which can be formulated in spatial, temporal, or spatiotemporal forms, are widely adopted to characterize non-Fickian diffusion and long-range particle transport in porous media [2, 8]. However, such single- or multi-term models are restricted to dynamics governed by a single or a finite set of fixed scaling exponents. These models produce mean-squared displacement following pure power-law forms and are incapable of reproducing logarithmic time-dependent diffusion behaviors. In complex multifractal media, fixed scaling exponents are rarely observed. Instead, the mean-squared displacement (MSD) often exhibits multiscale features or logarithmic time dependence. Consequently, distributed-order differential equations have been introduced as a powerful framework to model anomalous diffusion featuring a continuous superposition of scaling exponents [9, 10].

Differential equations incorporating distributed-order derivatives were first introduced by [11] to generalize stress–strain relationships in inelastic media. Subsequently, physical modeling paradigms for spatiotemporal distributed-order equations were systematically elaborated, highlighting that temporal distributed-order derivatives effectively capture memory effects, whereas spatial counterparts naturally account for media heterogeneity and non-local fluid-medium interactions [12–14]. Building upon these foundations, distributed-order kinetic theory was established by [10], revealing the statistical mechanisms of distributed-order operators in describing accelerating or decelerating subdiffusion and superdiffusion. In recent years, distributed-order models have been widely applied to wave propagation in viscoelastic rods [15, 16], ultraslow diffusion [17, 18, 19], Langevin stochastic dynamics [20], dynamic system modeling [21], and non-equilibrium transport [22]. Since closed-form analytical solutions are rarely accessible due to complex continuous integrals involving special functions, the development of robust and efficient numerical algorithms has become a central research focus. Explicit strong solutions and stochastic analogs on bounded domains were investigated by [9], while fundamental solutions and well-posedness were established via integral transforms by [17] and [22]. Regarding numerical schemes, general discretization frameworks were proposed by [23] and [24]. Implicit and compact difference schemes for temporal distributed-order equations were developed by [25], [26], [27], and [28]. Furthermore, spatial distributed-order models were studied on irregular domains by [29], while finite volume methods were formulated by [30, 31]. Additionally, second-order implicit methods based on shifted Grünwald–Letnikov operators and Crank–Nicolson strategies were developed by [1] and [2], respectively.

In summary, most existing numerical schemes for Riesz space distributed-order advection–dispersion equations only achieve second-order accuracy in both time and space [1, 30]; the available high-order schemes only cover the pure diffusion case and remain second-order accurate in time. To date, there is no numerical method with fourth-order accuracy in both time and space for such equations that includes fractional advection terms of order α ∈ ( 0,1 ) and fractional dispersion terms of order β ∈ ( 1,2 ). High-order spatio-temporal schemes can achieve prescribed accuracy with coarser grids, effectively reducing the computational cost of large-scale long-time simulations, and thus are of great research significance and practical value.

Motivated by the above research status, this article constructs a finite difference scheme that is fourth-order accurate in both time and space for the one-dimensional Riesz space distributed-order advection–dispersion equation. The main work is summarized as follows:

• The distributed-order integrals are discretized via the midpoint quadrature rule, converting the original equation into a multi-term Riesz space fractional advection–dispersion equation.
• The fourth-order difference formula derived from the fractional central difference operator combined with the weighted and shifted approach is employed to approximate Riesz fractional derivatives of all orders, yielding a semi-discrete system with fourth-order spatial accuracy.
• The (2,2) Padé approximation is adopted for temporal discretization, resulting in a fully implicit scheme with fourth-order accuracy in both time and space.
• Combining matrix spectral analysis with the energy method, we rigorously prove the unconditional stability of the proposed scheme and derive its global convergence order of O ( τ 4 + h 4 + 1 / L 2 + 1 / J 2 ).

The remainder of this article is organized as follows. Section 2 presents the numerical method. The stability and convergence analysis of the proposed numerical scheme is provided in Section 3. Section 4 provides a series of numerical experiments. Section 5 concludes the whole work.

2 Numerical methods

2.1 Preliminary

The fractional Sobolev space F n + γ ( R ) is defined bywhere υ ̂ ( w ) represents the Fourier transform of the function υ ( x ), and it is of the form υ ̂ ( w ) = ∫ R υ ( x ) e − i w x d x.

F n + γ R = υ | υ ∈ L 1 R , a n d ∫ R 1 + | w | n + γ | υ ̂ w | d w < ∞ ,

[32] proposed a fractional central difference operator ( ξ > − 1 ):whose coefficients are given by g k ( ξ ) = ( − 1 ) k Γ ( ξ + 1 ) Γ ( ξ / 2 ) − k + 1 Γ ( ξ / 2 ) + k + 1 and satisfy the following properties: g 0 ( ξ ) ≥ 0 and g k ( ξ ) = g − k ( ξ ) for all | k | ≥ 1.

L 0 , h ξ u x = ∑ k = − ∞ ∞ g k ξ u x − k h ,

[2] further proved that it provides a second-order approximation to the Riesz fractional derivative ( 1 < ξ ≤ 2 ):

∂ ξ u x ∂ | x | ξ = − 1 h ξ L 0 , h ξ u x + O h 2 .

The spatiotemporal domain [ 0 , T ] × [ a , b ] is discretized using uniform grids: t n = n τ for 0 ≤ n ≤ N and x i = a + i h for 0 ≤ i ≤ M, where τ = T / N and h = ( b − a ) / M. Based on the weighted and shifted difference idea, [33] constructed a fourth-order approximation formula for the Riesz fractional derivative, which is presented in the following lemma.

Lemma 2.1 [33] Suppose that u ( x ) is defined on [ a , b ] and belongs to F 4 + ξ ( R ) ξ ∈ ( 0,1 ) ⋃ ( 1,2 ) after zero extension to ( − ∞ , a ) ∪ ( b , + ∞ ) . Then, Moreover,

2.2 Numerical discretization of the Riesz space distributed-order advection–dispersion equation

Suppose that, for fixed ( x , t ), the mappingsbelong to C 2 [ 0,1 ] and C 2 [ 1,2 ], respectively. This regularity condition guarantees the second-order accuracy of the midpoint quadrature rule for the distributed-order integrals. We employ the midpoint quadrature rule [1] to discretize the distributed-order Riesz fractional derivative, which yieldswhere α l = 2 l − 1 2 L, l = 1,2 , … , L, and β j = 1 + 2 j − 1 2 J, j = 1,2 , … , J. Consequently, Problem (1) can be approximately converted into a multi-term Riesz fractional advection–dispersion equation.

α ↦ p α ∂ α u x , t ∂ | x | α , β ↦ q β ∂ β u x , t ∂ | x | β

∫ 0 1 p ⁡ α ∂ α u x , t ∂ | x | α ⁢ d α + ∫ 1 2 q ⁡ β ∂ β u x , t ∂ | x | β ⁢ d β = ∑ l = 1 L p α l L ⁢ ∂ α l u x , t ∂ | x | α l + ∑ j = 1 J q β j J ⁢ ∂ β j u x , t ∂ | x | β j + O ⁡ 1 L 2 + 1 J 2 ,

For Problem (Equation 1), assume that u ( ⋅ , t ) is defined on [ a , b ] and belongs to F 4 + β J ( R ) after zero extension to ( − ∞ , a ) ∪ ( b , + ∞ ). At the grid point ( x i , t ), applying the fourth-order approximation formula provided in Lemma 2.1 yieldswhere R i = O h 4 + 1 L 2 + 1 J 2. Denote u i ( t ) = u ( x i , t ) and f i ( t ) = f ( x i , t ), and define the vectors u ( t ) = [ u 1 ( t ) , … , u M − 1 ( t ) ] T, f ( t ) = [ f 1 ( t ) , … , f M − 1 ( t ) ] T, and R = [ R 1 , … , R M − 1 ] T. Then, Equation (2) can be rewritten in matrix-vector form aswith P = ∑ l = 1 L p ( α l ) L B ( α l ) h α l + ∑ j = 1 J q ( β j ) J B ( β j ) h β j, wherefor η ∈ { α l , β j } and all indices l , j. The matrices G − 1 η, G 0 η, and G 1 η are given as follows:

∂ u x i , t ∂ t = ∑ l = 1 L p α l L B x α l + ∑ j = 1 J q β j J B x β j u x i , t + f x i , t + R i , (2)

d u t d t = P u t + f t + R , (3)

B η = η 24 G − 1 η − 1 + η 12 G 0 η + η 24 G 1 η , (4)

G − 1 η = g 1 η g 0 η g − 1 η ⋯ g 5 − M η g 4 − M η g 3 − M η g 2 η g 1 η g 0 η g − 1 η ⋯ g 5 − M η g 4 − M η g 3 η g 2 η g 1 η g 0 η g − 1 η ⋯ g 5 − M η ⋮ ⋱ ⋱ ⋱ ⋱ ⋱ ⋮ g M − 3 η ⋯ g 3 η g 2 η g 1 η g 0 η g − 1 η g M − 2 η g M − 3 η ⋯ g 3 η g 2 η g 1 η g 0 η g M − 1 η g M − 2 η g M − 3 η ⋯ g 3 η g 2 η g 1 η ,

G 0 η = g 0 η g − 1 η g − 2 η ⋯ g 4 − M η g 3 − M η g 2 − M η g 1 η g 0 η g − 1 η g − 2 η ⋯ g 4 − M η g 3 − M η g 2 η g 1 η g 0 η g − 1 η g − 2 η ⋯ g 4 − M η ⋮ ⋱ ⋱ ⋱ ⋱ ⋱ ⋮ g M − 4 η ⋯ g 2 η g 1 η g 0 η g − 1 η g − 2 η g M − 3 η g M − 4 η ⋯ g 2 η g 1 η g 0 η g − 1 η g M − 2 η g M − 3 η g M − 4 η ⋯ g 2 η g 1 η g 0 η ,

G 1 η = g − 1 η g − 2 η g − 3 η ⋯ g 3 − M η g 2 − M η g 1 − M η g 0 η g − 1 η g − 2 η g − 3 η ⋯ g 3 − M η g 2 − M η g 1 η g 0 η g − 1 η g − 2 η g − 3 η ⋯ g 3 − M η ⋮ ⋱ ⋱ ⋱ ⋱ ⋱ ⋮ g M − 5 η ⋯ g 1 η g 0 η g − 1 η g − 2 η g − 3 η g M − 4 η g M − 5 η ⋯ g 1 η g 0 η g − 1 η g − 2 η g M − 3 η g M − 4 η g M − 5 η ⋯ g 1 η g 0 η g − 1 η .

By the properties of the coefficients g k ( η ), the matrix B ( η ) in Equation 4 is symmetric, so the matrix P in Equation 3 is also symmetric.

Suppose that the exact solution satisfies u ( x , ⋅ ) ∈ C 5 [ 0 , T ] with respect to t, which ensures that the temporal truncation error of the (2, 2) Padé approximation is of order O ( τ 5 ). Applying the (2,2) Padé approximation [34] to Equation 3 at time level t n yieldswhere u n = u ( t n ), f n = f ( t n ), Q n = d f n d t, and R ̂ i n = O τ ( τ 4 + h 4 + 1 L 2 + 1 J 2 ).

12 I − 6 τ P + τ 2 P 2 u n + 1 = 12 I + 6 τ P + τ 2 P 2 u n + 6 τ I − τ 2 P f n + 1 + 6 τ I + τ 2 P f n + τ 2 Q n − Q n + 1 + R ̂ n ,

Omitting the truncation error term R ̂ n, we arrive at the fully discrete numerical schemewhere W = ( 12 I − 6 τ P + τ 2 P 2 ) − 1 ( 12 I + 6 τ P + τ 2 P 2 ).

U n + 1 = W U n + 12 I − 6 τ P + τ 2 P 2 − 1 6 τ I − τ 2 P f n + 1 + 6 τ I + τ 2 P f n + τ 2 Q n − Q n + 1 ,

3 Stability and convergence

To derive the stability and convergence results for the proposed numerical scheme, we state several essential lemmas.

Lemma 3.1 [35] Let T n be a Toeplitz matrix with generating function, where f ( x ) is a 2 π -periodic continuous real-valued function. Denote λ min ( T n ) and λ max ( T n ) as the smallest and largest eigenvalues of T n , respectively. Then, we have where f min ( x ) and f max ( x ) denote the minimum and maximum values of f ( x ), respectively. Moreover, if f min ( x ) ≠ f max ( x ), then all eigenvalues of T n satisfy

Lemma 3.2 [33] Let β ∈ ( 1,2 ) ; the matrix − B ( β ) given in Equation 4 is real, symmetric, and positive definite.

Lemma 3.3 Let α ∈ ( 0,1 ) ; the matrix − B ( α ) given in Equation 4 is real, symmetric, and positive definite. Proof. The proof is completely analogous to Theorem 8 in [33]. The generating functions of the matrices G − 1 ( α ), G 0 ( α ), and G 1 ( α ) are given byUsing the identity ∑ k = − ∞ ∞ g k ( α ) e ikx = 2 ⁡ sin x 2 α, we obtain the generating function of the matrix − B ( α ) as follows One can readily observe that the function f − B ( α ) ( x ) is non-negative for ( α , x ) ∈ ( 0,1 ) × [ − π , π ]. By Lemma 3.1, the matrix − B ( α ) is positive definite, and hence, it is real, symmetric, and positive definite.

Lemma 3.4 The matrix P given in Equation 3 is real, symmetric, and negative definite. Proof. This result follows directly from the definition of matrix P together with Lemmas 3.2–3.3.

Lemma 3.5 For the matrix P in Equation 3, the following inequalities hold:where ‖ ⋅ ‖ 2 stands for the spectral norm. Proof. Take an arbitrary vector v ∈ R M − 1. Since P is real, symmetric, and negative definite, we have Substitute ( 12 I − 6 τ P + τ 2 P 2 ) − 1 v for v in inequality Equation 6. Noting that 12 I + 6 τ P + τ 2 P 2 commutes with its inverse, we obtainwhich yields ‖ ( 12 I − 6 τ P + τ 2 P 2 ) − 1 ( 12 I + 6 τ P + τ 2 P 2 ) ‖ 2 ≤ 1. Similarly, for any vector v ∈ R M − 1, we obtain Replacing v by ( 12 I − 6 τ P + τ 2 P 2 ) − 1 v in Equation 7, it follows thatwhich implies ‖ ( 12 I − 6 τ P + τ 2 P 2 ) − 1 ‖ 2 ≤ 1 12. This completes the proof.

Let V h = { v ∣ v = { v i } is a grid function defined on { x i } 0 M with v 0 = v M = 0 }. For v ∈ V h, the corresponding discrete L 2-norm is defined as

‖ v ‖ L 2 = h ∑ i = 1 M − 1 | v i | 2 1 2 .

Theorem 3.1 The numerical scheme Equation 5 is unconditionally stable.

Proof. Consider the perturbation equation. For an initial perturbation ε 0, its evolution satisfiesApplying the recurrence relation (Equation 8) n times, we obtain

ε n + 1 = W ε n , 0 ≤ n ≤ N − 1 .

ε n + 1 = W n + 1 ε 0 . (9)

Applying the vector Euclidean norm ‖ ⋅ ‖ to both sides of Equation 9 and using the submultiplicative property of matrix spectral norm yields

‖ ε n + 1 ‖ = ‖ W n + 1 ε 0 ‖ ≤ ‖ W ‖ 2 n + 1 ‖ ε 0 ‖ ≤ ‖ ε 0 ‖ .

Theorem 3.2 Suppose the following. belong to C 2 [ 0,1 ] and C 2 [ 1,2 ], respectively. Then, the numerical scheme Equation 5 is convergent with order O ( τ 4 + h 4 + 1 L 2 + 1 J 2 ).

By the recurrence relation and noting that the initial error satisfies e 0 = 0, we obtain

e n + 1 = ∑ i = 0 n W i K R ̂ n − i . (11)

Applying the L 2-norm to both sides of Equation 11 and using the submultiplicative property of matrix spectral norm yields

‖ e n + 1 ‖ L 2 = ‖ ∑ i = 0 n W i K R ̂ n − i ‖ L 2 ≤ ‖ K ‖ 2 ∑ i = 0 n ‖ W ‖ 2 i ‖ R ̂ n − i ‖ L 2 .

From Lemma 3.5, we have ‖ K ‖ 2 ≤ 1 12 and ‖ W ‖ 2 ≤ 1. It follows thatCombined with the local truncation error estimate R ̂ i n = O τ ( τ 4 + h 4 + 1 L 2 + 1 J 2 ) and noting that 0 ≤ ( n + 1 ) τ ≤ T, the total number of time steps satisfies n + 1 = O ( 1 τ ), which cancels the factor τ in the local truncation error. Hence, there exists a constant C independent of τ, h, 1 L, and 1 J, such that

‖ e n + 1 ‖ L 2 ≤ 1 12 ∑ i = 0 n ‖ R ̂ n − i ‖ L 2 .

‖ e n + 1 ‖ L 2 ≤ C τ 4 + h 4 + 1 L 2 + 1 J 2 .

4 Numerical experiments

In this section, we design a set of numerical experiments to verify the convergence orders of the proposed numerical scheme with respect to the distributed-order error, spatial error, and temporal error. Comparisons are also made with the existing finite difference method [1] to demonstrate the advantages of the proposed scheme.

Example 4.1 Consider the following Riesz space distributed-order advection–dispersion equation The exact solution is given by u ( x , t ) = e t x m ( 1 − x ) m with m ≥ 4, m ∈ N, whereand the source term readswhere in which c m , k , i denote the coefficients of the polynomial ∏ s = m + k + 1 2 m ( s − ξ ) = ∑ i = 0 m c m , k , i ξ i.

Set m = 4. We fix the spatial and temporal step sizes as h = τ = 1 / 100. By varying the step sizes of the midpoint quadrature rule, we compute the numerical errors and corresponding convergence orders. The results are reported in Table 1, which demonstrates that the midpoint quadrature rule achieves second-order accuracy.

TABLE 1. Distributed-order error results for Example 4.1 at t = 1 ( h = τ = 1 100 ).

Subsequently, we fix the step sizes of the midpoint quadrature rule such that 1 / L = 1 / J = 1 / 4000 and set the time step τ = 1 / 200. Numerical errors are computed by refining the spatial step size, and the resulting convergence orders are presented in Table 2. It can be observed that the spatial convergence order is four.

TABLE 2. Spatial errors and spatial convergence orders for Example 4.1 at t = 1 with different m 1 L = 1 J = 1 4000 and τ = 1 200.

Next, we keep 1 / L = 1 / J = 1 / 4000 for the midpoint quadrature rule and impose h = τ. We refine the spatial step size to obtain numerical errors and measure the associated convergence orders, as summarized in Table 3. The numerical results indicate that the temporal convergence order also attains fourth-order accuracy.

TABLE 3. Numerical errors and temporal convergence orders for Example 4.1 at t = 1 1 L = 1 J = 1 4000 and h = τ.

Finally, we compare our proposed scheme with the finite difference scheme given in [1]. The comparative results are listed in Table 4. Under identical step-size configurations, our scheme yields smaller numerical errors than the scheme from [1]. When the scheme in [1] is implemented on finer grids to reduce its errors to a magnitude comparable to that of our method, the CPU time increases drastically. Collectively, these numerical experiments verify the effectiveness of our proposed approach.

TABLE 4. Numerical errors and temporal convergence orders for Example 4.1 at t = 1 1 L = 1 J = 1 4000 and h = τ.

Remark 4.1 In the convergence theorem, the condition u ( ⋅ , t ) ∈ F 4 + β J ( R ) is sufficient for obtaining a fourth-order spatial error. We have tested several cases with different levels of spatial regularity, and the numerical results are presented in Table 2 . It can be observed from the table that when m ≥ 4 , the fourth-order convergence rate is indeed achieved; however, when the regularity requirement is lower (i.e., for m < 4 ), the theoretical fourth-order accuracy is not attained.

Remark 4.2 (Computational complexity and storage). We first analyze the cost of the direct implementation used in our numerical experiments. The matrix P is a weighted sum of L + J Toeplitz matrices. Assembling P explicitly costs O ( ( L + J ) N x 2 ) , where N x = M − 1 is the number of interior grid points. Computing P 2 , assembling A = 12 I − 6 τ P + τ 2 P 2 , and performing the LU factorization all require O ( N x 3 ) operations. Since A is time-independent, the LU decomposition is performed only once. Each subsequent time step involves forward and backward substitutions at a cost of O ( N x 2 ) , leading to a total complexity of O ( ( L + J ) N x 2 + N x 3 + N t N x 2 ) and a memory requirement of O ( N x 2 ). We emphasize that the above cost corresponds to a naive dense solver. Since B ( η ) are Toeplitz matrices, one may instead adopt a preconditioned conjugate gradient method with FFT-accelerated matrix-vector products (zero-padding is commonly used to suppress circulant wrap-around artifacts for Dirichlet-truncated Toeplitz matrices), where each B ( η ) v costs O ( N x ⁡ log N x ) . A suitable sine-transform-based preconditioner is expected to keep the iteration count nearly independent of N x , reducing the per-step cost to O ( ( L + J ) N x ⁡ log N x ) and the storage to O ( ( L + J ) N x ). This fast version eliminates the O ( N x 3 ) LU factorization overhead and is particularly attractive for large-scale simulations. A detailed investigation of this acceleration strategy is left for future work.

5 Conclusion

In this article, we propose an efficient numerical scheme with fourth-order accuracy in both time and space for the Riesz space distributed-order advection–dispersion equation. The midpoint quadrature rule is adopted to discretize the Riesz distributed-order derivative, which yields a multi-term Riesz space fractional advection–dispersion equation. A weighted and shifted fourth-order finite difference scheme is employed to discretize the Riesz space fractional advection–dispersion terms, resulting in the semi-discrete scheme. Afterward, the (2, 2) Padé approximation is utilized to discretize the temporal partial derivative, and thus, the fully discrete numerical scheme is obtained. The stability and convergence of the numerical scheme are rigorously established by combining spectral analysis with the energy method. Finally, a series of numerical experiments are designed to verify the effectiveness of the proposed numerical scheme.

Generative AI statement

The author(s) declared that generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Footnotes

Edited by: Jisheng Kou, Shaoxing University, China

Reviewed by: Xiangyi Peng, Xiangtan University, China

Fugui Ma, Peking University, China