Research article Special Issues

Efficient numerical schemes for variable-order mobile-immobile advection-dispersion equation

  • In this work, we present a high-order discontinuous Galerkin (DG) method with generalized alternating numerical fluxes to solve the variable-order (VO) fractional mobile-immobile advection-dispersion equation. This equation models complex transport phenomena where the order of differentiation varies with time, providing a more accurate representation of anomalous diffusion in heterogeneous media. For spatial and temporal discretization, the method employs the DG scheme and a finite difference method, respectively. Rigorous analysis confirms that the numerical scheme is unconditionally stable and convergent. Finally, numerical experiments are conducted to validate the theoretical results and illustrate the accuracy and efficiency of the scheme.

    Citation: Leqiang Zou, Yanzi Zhang. Efficient numerical schemes for variable-order mobile-immobile advection-dispersion equation[J]. Networks and Heterogeneous Media, 2025, 20(2): 387-405. doi: 10.3934/nhm.2025018

    Related Papers:

    [1] Raimund Bürger, Stefan Diehl, M. Carmen Martí, Yolanda Vásquez . A difference scheme for a triangular system of conservation laws with discontinuous flux modeling three-phase flows. Networks and Heterogeneous Media, 2023, 18(1): 140-190. doi: 10.3934/nhm.2023006
    [2] Changli Yuan, Mojdeh Delshad, Mary F. Wheeler . Modeling multiphase non-Newtonian polymer flow in IPARS parallel framework. Networks and Heterogeneous Media, 2010, 5(3): 583-602. doi: 10.3934/nhm.2010.5.583
    [3] Jinhu Zhao . Natural convection flow and heat transfer of generalized Maxwell fluid with distributed order time fractional derivatives embedded in the porous medium. Networks and Heterogeneous Media, 2024, 19(2): 753-770. doi: 10.3934/nhm.2024034
    [4] Zhe Pu, Maohua Ran, Hong Luo . Effective difference methods for solving the variable coefficient fourth-order fractional sub-diffusion equations. Networks and Heterogeneous Media, 2023, 18(1): 291-309. doi: 10.3934/nhm.2023011
    [5] Kexin Li, Hu Chen, Shusen Xie . Error estimate of L1-ADI scheme for two-dimensional multi-term time fractional diffusion equation. Networks and Heterogeneous Media, 2023, 18(4): 1454-1470. doi: 10.3934/nhm.2023064
    [6] Ming-Kai Wang, Cheng Wang, Jun-Feng Yin . A second-order ADI method for pricing options under fractional regime-switching models. Networks and Heterogeneous Media, 2023, 18(2): 647-663. doi: 10.3934/nhm.2023028
    [7] Patrick Henning, Mario Ohlberger . The heterogeneous multiscale finite element method for advection-diffusion problems with rapidly oscillating coefficients and large expected drift. Networks and Heterogeneous Media, 2010, 5(4): 711-744. doi: 10.3934/nhm.2010.5.711
    [8] Giuseppe Maria Coclite, Lorenzo di Ruvo . A singular limit problem for conservation laws related to the Kawahara-Korteweg-de Vries equation. Networks and Heterogeneous Media, 2016, 11(2): 281-300. doi: 10.3934/nhm.2016.11.281
    [9] Jin Cui, Yayun Fu . A high-order linearly implicit energy-preserving Partitioned Runge-Kutta scheme for a class of nonlinear dispersive equations. Networks and Heterogeneous Media, 2023, 18(1): 399-411. doi: 10.3934/nhm.2023016
    [10] Diandian Huang, Xin Huang, Tingting Qin, Yongtao Zhou . A transformed $ L1 $ Legendre-Galerkin spectral method for time fractional Fokker-Planck equations. Networks and Heterogeneous Media, 2023, 18(2): 799-812. doi: 10.3934/nhm.2023034
  • In this work, we present a high-order discontinuous Galerkin (DG) method with generalized alternating numerical fluxes to solve the variable-order (VO) fractional mobile-immobile advection-dispersion equation. This equation models complex transport phenomena where the order of differentiation varies with time, providing a more accurate representation of anomalous diffusion in heterogeneous media. For spatial and temporal discretization, the method employs the DG scheme and a finite difference method, respectively. Rigorous analysis confirms that the numerical scheme is unconditionally stable and convergent. Finally, numerical experiments are conducted to validate the theoretical results and illustrate the accuracy and efficiency of the scheme.



    Fractional differential equations (FDEs) have become an essential mathematical framework for modeling a wide range of complex phenomena, including long-range memory effects, anomalous mechanical dynamics, and advanced control processes. By generalizing classical calculus, fractional calculus enables a more accurate and flexible representation of physical systems that are inadequately captured by traditional integer-order models. This mathematical approach has been successfully applied in various disciplines, including quantitative finance, engineering, biology, chemistry, and hydrology [1,2]. Among recent developments, variable-order (VO) fractional calculus has garnered significant attention due to its capacity to describe systems with time- or space-dependent memory characteristics and evolving dynamical behavior [3].

    In recent years, a wide variety of numerical methods have been developed to solve fractional partial differential equations (FPDEs), which are often analytically intractable due to the nonlocal nature of fractional operators. Among these, finite element methods have been widely applied and rigorously analyzed for various types of time- and space-fractional models [4,5,6,7,8]. Finite volume methods have been adopted to handle anomalous transport and advection-diffusion processes with fractional dynamics [9,10,11]. Finite difference methods, known for their simplicity and effectiveness, have also been extensively investigated, with significant progress in high-order, compact, and implicit schemes [12,13,14,15,16]. Spectral and spectral-Galerkin methods offer highly accurate solutions and have demonstrated spectral convergence for a variety of FPDEs [17,18,19]. Meshless methods such as radial basis function and local interpolation approaches have emerged as powerful alternatives for multidimensional and complex-geometry problems [20,21,22,23]. Approximate analytical methods provide efficient solutions by combining analytical and numerical approaches, often using perturbation or homotopy techniques to simplify complex fractional equations [24,25,26]. To improve flexibility and accuracy, discontinuous Galerkin methods have been developed for both constant- and variable-order fractional equations [27,28,29,30].

    The discontinuous Galerkin (DG) method combines the advantages of finite element and finite volume methods. By employing a discontinuous solution space, the DG method provides high accuracy and flexibility in solving a wide range of partial differential equations, including those with complex geometries or discontinuous solutions. This method is particularly effective for approximations of arbitrary order, making it a preferred choice for a variety of computational problems. A key strength of the DG method lies in its ability to handle adaptive meshing and high-order polynomial approximations. Furthermore, the method is supported by rigorous theoretical analysis, including error estimates and stability results.

    The time-variable-order (VO) fractional mobile-immobile advection-dispersion equation effectively captures solute transport processes in porous and fractured media, providing a robust framework for studying complex migration phenomena [31,32,33]. A comprehensive overview of the physical background and formulation of this model is provided in [34,35]. Sadri et al. [36] employed a spectral collocation method utilizing sixth-kind Chebyshev polynomials to solve the VO time-fractional mobile-immobile advection equations, demonstrating the efficiency and accuracy of the method. Ma et al. [37] introduced the Jacobi spectral collocation method for the variable-order advection-dispersion equation, validating its higher-order convergence through numerical analysis. Liu et al. [38] analyzed a second-order finite difference scheme for fractal mobile-immobile transport equations. Golbabai et al. [39] applied meshless methods with radial basis functions to approximate solutions for the time-fractional mobile-immobile advection-dispersion model within bounded regions. Zhang et al. [40] proposed an implicit Euler method to solve the mobile-immobile advection-dispersion model, rigorously demonstrating its unconditional stability. Jiang et al. [41] developed a numerical method based on reproducing kernel theory combined with the collocation method to address the mobile-immobile advection-dispersion model. Saffarian and Mohebbi [42] developed a robust numerical scheme for the two-dimensional time-variable-order fractional mobile-immobile advection-dispersion model, proving the unconditional stability of the fully discrete scheme.

    In this paper, we propose and analyze a DG method based on generalized alternating numerical fluxes to solve the following variable-order mobile-immobile advection-dispersion equaiton

    c1ωt+c2Pα(t)tω=d1ωx+d2ωxx+g(x,t),(x,t)(a,b)×(0,T],ω(x,0)=ω0(x),x[a,b], (1.1)

    where 0<ρ(t)<1, c1>0, c2>0, and d1,d2>0. In this paper, g and ω0 are assumed to be smooth functions, and the solution is considered to be either periodic or compactly supported.

    The variable-order (VO) fractional derivative operator [43] is defined in the Caputo sense

    Pρ(t)tω(x,t)=1Γ(1ρ(t))t0ω(x,s)s1(ts)ρ(t)ds,

    where 0<ρ(t)<1.

    The mobile-immobile advection-dispersion equation is a fundamental tool for describing solute transport, yet its standard formulation assumes constant-order derivatives, which may not accurately reflect the varying memory effects and heterogeneous diffusion characteristics observed in real-world systems. To address these limitations, we employ a variable-order (VO) fractional derivative, allowing the diffusion process to dynamically adapt to spatial and temporal variations in the medium. The main contributions of this paper are as follows:

    (1) A high-order DG scheme with generalized alternating numerical fluxes is presented to efficiently solve the VO fractional mobile-immobile advection-dispersion equation.

    (2) The stability and convergence of the proposed scheme are analyzed, ensuring its reliability for practical applications.

    (3) Extensive numerical experiments are conducted to validate the theoretical findings and demonstrate the accuracy of the method for the variable-order equation.

    In Section 2, the symbols, projections, and theorems required for the proof are introduced. Section 3 presents the construction of the numerical scheme for Eq (1.1) using the LDG method, and it is demonstrated that the scheme is unconditionally stable and convergent with an accuracy of O(Δt+hk+1). In Section 4, numerical experiments are conducted to validate the effectiveness and reliability of the proposed scheme. Finally, the conclusions are summarized in Section 5.

    The interval Ω=[a,b] is partitioned into subintervals, denoted by J, where a=x12<x32<<xN+12=b. For j=1,,N, each subinterval is defined as Ij=[xj12,xj+12], with the element size given by Δxj=xj+12xj12. The maximum element size is denoted by h=max1jNΔxj.

    The time interval [0,T] is discretized uniformly into M time steps, such that Δt=TM. The discrete time points are given as tn=nΔt for n=0,1,,M.

    At each boundary point xj+12, the left and right limits of the function ω are defined as ωj+12 and ω+j+12, respectively. Specifically, ωj+12 corresponds to the value in the left cell Ij, while ω+j+12 corresponds to the value in the right cell Ij+1.

    The discontinuous Galerkin space Vk is defined as:

    Vk={vPk(Ij):xIj,j=1,2,,N},

    where Pk(Ij) denotes the space of polynomials of degree k defined on each subinterval Ij.

    In the paper, C is used to represent a positive constant, which may take different values in different cases. The inner product on L2(D) is denoted by (,)D, and the associated norm is represented by D. For D=Ω, the subscript is omitted for simplicity.

    First, we introduce some lemmas that will be used to design the scheme.

    Lemma 3.1. [1] For 0<ρ(t)<1 and t>0, the left-sided Caputo fractional derivative is defined as

    C0Pρ(t)th(x,t)=1Γ(1ρ(t))t0h(x,s)s1(ts)ρ(t)ds,

    is equivalent to the Riemann-Liouville fractional derivative, expressed as

    Pρ(t)th(x,t)=1Γ(1ρ(t))ddtt0h(x,s)(ts)ρ(t)ds,

    provided that the condition h(x,0)=0 is satisfied.

    Lemma 3.2. [1] (Grünwald formula)

    The function space Ψm+a(R) is defined as

    Ψm+a(R)={ω | ωL1(R):(1+|ξ|)m+a|ˆω(ξ)|dξ<}.

    Let ωΨ1+a(R). It follows that

    1Γ(1ρ(t))ddttω(s)(ts)ρ(t)ds=1(Δt)ρ(t)j=0ρρ(t)jω(t(jr)Δt)+O(Δt),

    where r is an integer, and θρ(t)j is given by

    θρ(t)j=(1)j(ρ(t)j).

    The coefficients θρ(t)j satisfy the following properties for 0<ρ(t)<1:

    θρ(t)0=1,θρ(t)1=ρ(t)0,
    θρ(t)2θρ(t)3θρ(t)40,
    k=0θρ(t)k=0,nk=0θρ(t)k0for n1.

    Moreover, these coefficients can be evaluated recursively as

    θρ(t)0=1,θρ(t)k=(1ρ(t)+1k)θρ(t)k1for k1.

    Lemma 3.3. There is a constant σ>0 such that the following condition holds:

    c2(Δt)ρ(t)n1k=0θρ(t)kσ>0.

    Let M be a positive integer, and define tn=nMT. The values of ωt and Pρ(tn)tu at tn are approximated as follows:

    ωt(x,tn)=ω(x,tn)ω(x,tn1)Δt+γn1(x), (3.1)

    where γn1(x) represents the truncation error in the approximation of ωt. Additionally, by applying Lemmas 3.1 and 3.2, it is obtained that

    Pρ(tn)tω(x,tn)=1Γ(1ρ(tn))tn0ω(x,s)s1(tns)ρ(tn)ds=1(Δt)ρ(tn)nk=0θρ(tn)kω(x,(nk)Δt)+γn2(x), (3.2)

    where γn2(x) denotes the truncation error associated with the approximation of Pρ(tn)tω. The total truncation error in the temporal direction is given by γn(x)=γn1(x)+γn2(x), and it satisfies the bound γn(x)CΔt.

    The problem (1.1) can be reformulated into the following system:

    p=ωx,c1ωt+c2Pρ(t)tω+d1ωxd2px=g(x,t). (3.3)

    Let ωnh,pnhVk denote the approximations of ω(,tn) and p(,tn), respectively, and define gn(x)=g(x,tn). The fully discrete scheme is then formulated as follows: find ωnh,pnhVk such that, for all v,μVk,

    c1ΔtΩωnhvdx+c2(Δt)ρ(tn)θρ(tn)0Ωωnhvdx+d1(Ωωnhvxdx+Nj=1((~ωnhv)j+12(~ωnhv+)j12))+d2(ΩpnhvxdxNj=1((^pnhv)j+12(^pnhv+)j12))=c1ΔtΩωn1hvdx+c2(Δt)ρ(tn)nk=1(θρ(tn)k)Ωωnkhvdx+Ωgnvdx,Ωpnhμdx+ΩωnhμxdxNj=1((^ωnhμ)j+12(^ωnhμ+)j12)=0. (3.4)

    The hat functions in the boundary terms, which arise from the integration by parts in Eq (3.4), are referred to as numerical fluxes. The choice of an appropriate numerical flux plays a critical role in the theoretical analysis of the LDG scheme. From a practical perspective, generalized alternating numerical fluxes offer greater flexibility and broader applicability compared to traditional numerical fluxes [44]. The generalized alternating numerical fluxes are defined as follows:

    ^ωnh=δ(ωnh)+(1δ)(ωnh)+,^pnh=(1δ)(pnh)+δ(pnh)+,~ωnh=δ(ωnh)+(1δ)(ωnh)+, (3.5)

    where δ[0,12)(12,1]. For δ=12, the properties related to the uniqueness and approximation of the generalized Gauss-Radau projection become more intricate.

    For convenience, we denote

    ΦΩ(ωnh,pnh;μ,v)=ΩωnhμxdxNj=1(((ωnh)μ)j+12((ωnh)μ+)j12)+ΩpnhvxdxNj=1(((pnh)+v)j+12((pnh)+v+)j12). (3.6)

    To simplify the analysis, we focus on the case where g=0 in the theoretic analysis.

    Theorem 3.1. Suppose that ω(x,t)C([0,T],Hs(Ω)) with sk+1. Then, the fully discrete LDG scheme (3.4), with the flux (3.5), is unconditionally stable. Furthermore, ωnh satisfies

    ωnhω0h,n=1,2,,M. (3.7)

    Proof. Taking v=ωnh and μ=d2pnh in Eq (3.4) and using the flux choice in Eq (3.5), we derive the following equation:

    c1Δtωnh2+c2(Δt)ρ(tn)θρ(tn)0ωnh2+12d1Nj=1[ωnh]2j12+d2pnh2+ΦΩ(ωnh,pnh;d2pnh,ωnh)=c1ΔtΩωn1hωnhdx+c2(Δt)ρ(tn)nk=1(θρ(tn)k)Ωωnkhωnhdx. (3.8)

    Next, consider the cell Ij=[xj12,xj+12]. The term ΦIj(ωnh,pnh;d2pnh,ωnh) is computed as:

    ΦIj(ωnh,pnh;d2pnh,ωnh)=d2(Ijωnh(pnh)xdx((ωnh)(pnh))j+12+((ωnh)(pnh)+)j12)+d2(Ijpnh(ωnh)xdx((pnh)+(ωnh))j+12+((pnh)+(ωnh)+)j12)=d2(((pnh)(ωnh))j+12((pnh)+(ωnh)+)j12((ωnh)(pnh))j+12+((ωnh)(pnh)+)j12((pnh)+(ωnh))j+12+((pnh)+(ωnh)+)j12). (3.9)

    By summing over j=1,,N in Eq (3.9), we observe that all boundary terms cancel out. Therefore, we have:

    ΦΩ(ωnh,pnh;d2pnh,ωnh)=0. (3.10)

    By Lemma 3.2, we have

    θρ(tn)nn1k=0θρ(tn)k. (3.11)

    Combining Eq (3.10) and the Cauchy-Schwarz inequality, Eq (3.8) becomes

    c1Δtωnh2+c2(Δt)ρ(tn)θρ(tn)0ωnh2+12d1Nj=1[ωnh]2j12+d2pnh2c1Δtωn1hωnh+c2(Δt)ρ(tn)n1k=1(θρ(tn)k)ωnkhωnh+c2(Δt)ρ(tn)(θρ(tn)n)ω0hωnhc1Δtωn1hωnh+c2(Δt)ρ(tn)n1k=1(θρ(tn)k)ωnkhωnh+c2(Δt)ρ(tn)n1k=0θρ(tn)kω0hωnh, (3.12)

    where the last inequality uses Eq (3.11).

    Dividing both sides by ωnh, we obtain

    (c1Δt+c2(Δt)ρ(tn)θρ(tn)0)ωnhc1Δtωn1h+c2(Δt)ρ(tn)n1k=1(θρ(tn)k)ωnkh+c2(Δt)ρ(tn)n1k=0θρ(tn)kω0h. (3.13)

    To prove the theorem, we use mathematical induction.

    For n=1, Eq (3.13) simplifies to

    ω1hω0h. (3.14)

    Assume that

    ωmhω0h,m=1,2,,n1. (3.15)

    Using Eq (3.13), we need to prove that ωnhω0h. Substituting the inductive hypothesis into Eq (3.13), we have

    (c1Δt+c2(Δt)ρ(tn)θρ(tn)0)ωnhc1Δtω0h+c2(Δt)ρ(tn)n1k=1(θρ(tn)k)ω0h+c2(Δt)ρ(tn)n1k=0θρ(tn)kω0h=(c1Δt+c2(Δt)ρ(tn)θρ(tn)0)ω0h. (3.16)

    Dividing through by c1Δt+c2(Δt)ρ(tn)θρ(tn)0, we obtain

    ωnhω0h.

    By mathematical induction, the inequality holds for all n. Hence, the theorem is proved.

    First, we introduce the generalized Gauss–Radau projection, which will be used in proving convergence.

    For any periodic function ϖ defined on the interval [a,b], the generalized Gauss-Radau projection [44], denoted by Qδϖ, is uniquely determined. Let ϖe=Qδϖϖ be the corresponding projection error. When δ12, the following conditions are satisfied for j=1,2,,N:

    Ijϖevdx=0,vPk1(Ij),and(ϖe)(δ)j+12=0. (3.17)

    Based on these conditions, the following result can be established [44,45].

    Lemma 3.4. Let δ12. If ϖHs+1[a,b], the inequality below holds:

    ϖe+h12ϖeL2(Γh)Chmin(k+1,s+1)ϖs+1, (3.18)

    where C>0 is a constant independent of h and ϖ.

    Theorem 3.2. Let ω(x,t)C([0,T],Hs(Ω)), with sk+1, denote the exact solution of the problem (1.1), and let ωnh represent the numerical solution obtained using the LDG scheme (3.4). The following error estimate holds

    ω(x,tn)ωnhC(Δt+hk+1),

    where C is a constant that depends on u and T.

    Proof.

    enω=ω(x,tn)ωnh=ξnωηnω,ξnω=Qδenω,ηnω=Qδω(x,tn)ω(x,tn),enp=p(x,tn)pnh=ξnpηnp,ξnp=Q1δenp,ηnp=Q1δp(x,tn)p(x,tn). (3.19)

    The terms ηnω and ηnp are estimated using inequality (3.4).

    Starting from the flux expression (3.5), with the test functions v=ξnω and μ=d2ξnp, we obtain the following equation:

    c1ΔtΩ(ξnω)2dx+c2(Δt)ρ(tn)θρ(tn)0Ω(ξnω)2dx+d12Nj=1[ξnω]2j12+d2Ω(ξnp)2dx+ΦΩ(ξnω,ξnp;d2ξnp,ξnω)=c1ΔtΩηnωξnωdx+c2(Δt)ρ(tn)θρ(tn)0Ωηnωξnωdx+d2Ωηnpξnpdxc1ΔtΩηn1ωξnωdxΩγn(x)vdx+ΦΩ(ηnω,ηnp;d2ξnp,ξnω)+c2(Δt)ρ(tn)nk=1(θρ(tn)k)Ωξnkωξnωdxc2(Δt)ρ(tn)nk=1(θρ(tn)k)Ωηnkωξnωdx. (3.20)

    Using properties (3.17), it follows that

    ΦΩ(ηnω,ηnp;d2ξnp,ξnω)=0.

    Next, applying Eqs (3.10) and (3.11), and the Cauchy-Schwarz inequality, the following inequality is obtained

    c1Δtξnω2+c2(Δt)ρ(tn)θρ(tn)0ξnω2+d12Nj=1[ξnω]2j12+d2ξnp2c1Δt(ηnω+ηn1ω)ξnω+c1Δtξn1ωξnω+γn(x)ξnω+c2(Δt)ρ(tn)n1k=1(θρ(tn)k)ξnkωξnω+c2(Δt)ρ(tn)n1k=0θρ(tn)kξ0ωξnω+c2(Δt)ρ(tn)n1k=1(θρ(tn)k)ηnkωξnω+c2(Δt)ρ(tn)n1k=0θρ(tn)kη0ωξnω+c2(Δt)ρ(tn)θρ(tn)0ηnωξnω+d2ηnpξnp. (3.21)

    Using the estimate

    ηiωηi1ωΔt1Δttiti1t(Qδω(x,t)ω(x,t))Chk+1ωtL(H2(Ω)),

    we obtain

    Ω(ηnωηn1ωΔt)ξnωdxChk+1ωtL(H2(Ω))ξnω,

    and thus, we have the following inequality:

    c1Δtξnω2+c2(Δt)ρ(tn)θρ(tn)0ξnω2+d12Nj=1[ξnω]2j12+d2ξnp2(c1Δt(ηnω+ηn1ω)+c1Δtξn1ω+γn(x)+(c2(Δt)ρ(tn)n1k=1(θρ(tn)k))(ξnkω+ηnkω)+(c2(Δt)ρ(tn)n1k=0θρ(tn)k)(ξ0ω+η0ω)+c2(Δt)ρ(tn)θρ(tn)0ηnω)ξnω+d2ηnpξnp. (3.22)

    There exists a constant χ>0 such that

    (θρ(tn)1+θρ(tn)2++θρ(tn)n1)θρ(tn)0χ(θρ(tn)1+θρ(tn)2++θρ(tn)n1),

    using Lemma 3.3 and noticing the fact that

    abεa2+14εb2,aR,bR,ε>0. (3.23)

    Let k=c1Δt+c2(Δt)ρ(tn)θρ(tn)0, and then we use the fact (3.23) twice, first let

    a=(c1Δt(ηnω+ηn1ω)+c1Δtξn1ω+γn(x)+(c2(Δt)ρ(tn)n1k=1(θρ(tn)k))(ξnkω+ηnkω)+(c2(Δt)ρ(tn)n1k=0θρ(tn)k)(ξ0ω+η0ω)+c2(Δt)ρ(tn)θρ(tn)0ηnω),b=ξnω,ε=k2>0,

    then let a=ηnp,b=ξnp,ε=12. Since

    (c2(Δt)ρ(tn)n1k=1χ(θρ(tn)k))(c2(Δt)ρ(tn)θρ(tn)0),

    so

    c2(Δt)ρ(tn)θρ(tn)0ηnωχ(c2(Δt)ρ(tn)n1k=1(θρ(tn)k))ηnω,

    then the inequality (3.22) becomes

    k2ξnω212k(c1Δt(ηnω+ηn1ω)+c1Δtξn1ω+(c2(Δt)ρ(tn)n1k=1(θρ(tn)k))(ξnkω+ηnkω+χηnω)+(c2(Δt)ρ(tn)n1k=0θρ(tn)k)(ξ0ω+η0ω+1σγn(x)))2+d22ηnp2. (3.24)

    By the fact that

    a2+b2(a+b)2,ab0,

    multiplying Eq (3.24) by 2k, and let

    a=c1Δt(ηnω+ηn1ω)+c1Δtξn1ω+(c2(Δt)ρ(tn)n1k=1(θρ(tn)k))(ξnkω+ηnkω+χηnω)+(c2(Δt)ρ(tn)n1k=0θρ(tn)k)(ξ0ω+η0ω+1σγn(x)),b=d2kηnp,

    we can obtain the following inequality:

    kξnωc1Δt(ηnω+ηn1ω)+c1Δtξn1ω+(c2(Δt)ρ(tn)n1k=1(θρ(tn)k))(ξnkω+ηnkω+χηnω)+d2kηnp+(c2(Δt)ρ(tn)n1k=0θρ(tn)k)(ξ0ω+η0ω+1σγn(x)). (3.25)

    When n=1, Eq (3.25) becomes

    (c1Δt+c2(Δt)ρ(t1)θρ(t1)0)ξ1ωc1Δt(η1ω+η0ω)+c1Δtξ0ω+d2kη1p+(c2(Δt)ρ(t1)θρ(t1)0)(ξ0ω+η0ω+1σγ1(x)). (3.26)

    Note that ξ0ω=0, and from Eq (3.4), we obtain

    ξ1ωC(Δt+hk+1). (3.27)

    We assume that the inequality holds for all m=1,2,,n1, that is

    ξmωC(Δt+hk+1). (3.28)

    It is noted that

    c2(Δt)ρ(tn)n1k=0θρ(tn)k+c2(Δt)ρ(tn)n1k=1(θρ(tn)k)=c2(Δt)ρ(tn)θρ(tn)0.

    From Eqs (3.25) and (3.28), the following inequality is obtained:

    (c1Δt+c2(Δt)ρ(tn)θρ(tn)0)ξnωc1ΔtC(Δt+hk+1)+c2(Δt)ρ(tn)n1k=1(θρ(tn)k)C(Δt+hk+1)+c2(Δt)ρ(tn)n1k=0θρ(tn)kC(Δt+hk+1). (3.29)

    Thus, we conclude that

    ξnωC(Δt+hk+1).

    Finally, by applying the triangle inequality and the interpolation property, the proof of Theorem 3.2 is completed.

    Example 4.1. Consider the problem (1.1) with parameters c1=1, c2=2, d1=1, d2=2, and the initial condition u(x,0)=0. The function g(x,t) is defined as

    g(x,t)=2tsin(2πx)+4t2ρ(t)Γ(3ρ(t))sin(2πx)+2πt2cos(2πx)+4π2t2sin(2πx).

    The exact solution for this problem is given by ω(x,t)=t2sin(2πx).

    The spatial convergence properties of the scheme (3.4) are evaluated using this example. The temporal step size is fixed at Δt=1/1000, while the spatial mesh sizes are varied as h=1/5, 1/10, 1/15, and 1/20, respectively. The resulting numerical errors and convergence rates in both the L2-norm and L-norm, corresponding to various fractional orders, are presented in Tables 1 and 2. It is evident that the proposed scheme achieves the optimal convergence rates when using piecewise Pk polynomial.

    Table 1.  Spatial accuracy test when taking piecewise Pk polynomials using generalized numerical fluxes on uniform meshes for ρ(t)=et+15, Δt=1/1000,δ=0.6,T=1.
    ρ(t) Pk N L-error order L2-error order
    5 0.781356543644563 - 0.811643551435544 -
    10 0.390678271822282 1.00 0.431944597464271 0.91
    P0 15 0.260452181214854 0.92 0.298665460272027 0.92
    20 0.195339135911141 0.97 0.229874475006395 0.89
    5 0.365454345536453 - 0.413153543451345 -
    ρ(t)=et+15 10 0.093283329546276 1.97 0.109937078106140 1.91
    P1 15 0.041966645713413 1.92 0.050676881906131 1.94
    20 0.023810852648272 1.92 0.029253436970555 1.97
    5 0.025342434342468 - 0.032634552345326 -
    10 0.003466508284447 2.87 0.004341905554260 2.91
    P2 15 0.001082705465106 2.82 0.001334304063205 2.91
    20 0.000474172272630 2.91 0.000577674350873 2.92

     | Show Table
    DownLoad: CSV
    Table 2.  Spatial accuracy test when taking piecewise Pk polynomials using generalized numerical fluxes on uniform meshes for ρ(t)=sint+210, Δt=1/1000,δ=0.2,T=1.
    ρ(t) Pk N L-error order L2-error order
    5 0.543432507028769 - 0.732453465895323 -
    10 0.292479082134599 0.92 0.367236212865906 0.95
    P0 15 0.197158461423066 0.93 0.247942026683542 0.94
    20 0.148119129679174 0.96 0.187525639215507 0.90
    5 0.386745645464365 - 0.429856341235466 -
    ρ(t)=sint+210 10 0.092682012633183 1.87 0.099368444058785 1.98
    P1 15 0.039432065256851 1.82 0.045507256764390 1.94
    20 0.021259111655579 1.92 0.023036687278259 1.91
    5 0.035464563245658 - 0.042325479864524 -
    10 0.006003059928863 2.85 0.007044728419360 2.90
    P2 15 0.002187236794468 2.92 0.002564554191487 2.94
    20 0.001058292464021 2.91 0.001241791511304 2.92

     | Show Table
    DownLoad: CSV

    By fixing a sufficiently small spatial mesh size h=1500 and choosing various temporal step sizes Δt=15,110,120,140, we observe from Table 3 that the numerical results exhibit a first-order convergence in time, which is consistent with our theoretical results.

    Table 3.  Temporal errors and convergence rates for piecewise linear P1 basis functions when Δt=1/1000,δ=0.3,T=1.
    M L2-error order L-error order
    5 0.135345464564645 - 0.093453566456722 -
    10 0.067198735409812 1.01 0.046712374092112 1.03
    ρ(t)=(cost)2152t 20 0.032480178908534 1.04 0.022305479807943 1.05
    40 0.016432954718745 0.98 0.010942573201284 1.07
    5 0.134556576878564 - 0.114547748454358 -
    10 0.061872654328745 1.10 0.055431287635421 1.02
    ρ(t)=1+et10 20 0.029389751008234 1.05 0.026058974358762 0.97
    40 0.014032178945671 1.01 0.012462398754123 1.08

     | Show Table
    DownLoad: CSV

    To further validate the accuracy of the proposed high-order DG method, we perform a comparison with a finite difference (FD) method for solving the variable-order fractional advection-dispersion Eq (1.1) with the same parameters and the initial condition. We first divide the spatial domain [a,b] into N equal subintervals with mesh size h1=baN, and denote xj=a+jh1, for j=0,1,,N. We approximate the spatial derivatives of ω(x,tn) at the grid points xj using finite differences as follows:

    ωx(xj,tn)ω(xj+1,tn)ω(xj,tn)h1,ωxx(xj,tn)ω(xj+1,tn)2ω(xj,tn)+ω(xj1,tn)h21. (4.1)

    Substitute Eqs (4.1), (3.1) and (3.2) into the original equation, we obtain the fully discrete scheme

    c1Δt(ωnjωn1j)+c2(Δt)ρ(tn)nk=0θρ(tn)kωnkj=d1ωnj+1ωnj12h1+d2ωnj+12ωnj+ωnj1h21+gnj, (4.2)

    where

    ωnjω(xj,tn),

    gnj=g(xj,tn) is the known source term,

    θρ(tn)k are weights from the fractional convolution quadrature,

    γnj is the total truncation error at (xj,tn), satisfying γnjC(Δt+h21).

    The variable fractional order is chosen as ρ(t)=8t15.

    In the comparison, we fix the time step as Δt=1/1000 and vary the spatial mesh size h1=1/10,1/20,1/30,1/40. The numerical errors and observed convergence rates in the L2-norm for both the DG method (using piecewise P1 polynomials) and the FD method (4.2) are shown in Figure 1. It is observed that both methods exhibit second-order convergence in space now. However, the DG method demonstrates higher accuracy than the finite difference scheme. Moreover, the DG method is capable of achieving (k+1)-th order spatial convergence when using piecewise Pk polynomials. For example, third-order convergence is attained with piecewise P2 polynomials, as evidenced by the results in Tables 1 and 2.

    Figure 1.  Error in L2 norm for DG method (using piecewise P1 polynomials) and the FD method (4.2), ρ(t)=8t15.

    Example 4.2. To illustrate the impact of different variable-order fractional derivatives, we consider the initial condition:

    ω(x,0)=sin(πx),x[0,1].

    We solve the variable-order mobile-immobile advection-dispersion equation (1.1) using two different choices for the fractional order function:

    ● Case 1: A linear function, ρ(t)=0.80.3t.

    ● Case 2: An exponential function, ρ(t)=0.6e0.5t.

    The final time is set to T=1, and the computed solutions are plotted in Figure 2. We can observe

    ● The solution with linear order ρ(t)=0.80.3t leads to a gradual decrease in diffusion, resulting in a smoother and more spread-out profile.

    ● The solution with exponential order ρ(t)=0.6e0.5t exhibits initially strong diffusion that slows down over time, leading to a more localized concentration of the solution.

    ● This comparison highlights how different variable-order fractional derivatives influence the dispersion behavior, with lower orders leading to slower diffusion and more localized effects.

    Figure 2.  Numerical solutions at T=1 for different fractional orders.

    These results demonstrate the flexibility of variable-order fractional models in capturing complex transport phenomena, making them more suitable for describing real-world diffusion processes compared to classical integer-order models.

    In this study, a high-order numerical method is introduced for solving the time-variable-order fractional mobile-immobile advection-dispersion model. The proposed scheme combines the finite difference method with the local discontinuous Galerkin (LDG) method. Through the careful selection of projections and numerical fluxes, it is demonstrated that the scheme is unconditionally stable and achieves a convergence rate of O(Δt+hk+1) in the L2 norm.

    It is worth noting that while this work focuses on the time-fractional model, the method could also be extended to handle fractional derivatives with respect to the spatial variable, as explored in [46]. The numerical scheme and stability analysis would require modifications to accommodate the nonlocal nature of spatial fractional derivatives, which will be considered in our future work.

    L. Zou wrote the main manuscript text. Y. Zhang analyzed the stability and convergence of the scheme. All authors reviewed the manuscript.

    The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

    The authors declare there is no conflict of interest.



    [1] I. Podlubny, Fractional Differential Equations, Academic Press, 1999.
    [2] W. Deng, Smoothness and stability of the solutions for nonlinear fractional differential equations, Nonlinear Anal. Theory Methods Appl., 72 (2010), 1768–1777. https://doi.org/10.1016/j.na.2009.09.018 doi: 10.1016/j.na.2009.09.018
    [3] Q. Li, Y. Chen, Y. Huang, Y. Wang, Two-grid methods for nonlinear time fractional diffusion equations by L1-Galerkin FEM, Math. Comput. Simul., 185 (2021), 436–451. https://doi.org/10.1016/j.matcom.2020.12.033 doi: 10.1016/j.matcom.2020.12.033
    [4] X. Zheng, H. Wang, Optimal-order error estimates of finite element approximations to variable-order time-fractional diffusion equations without regularity assumptions of the true solutions, IMA J. Numer. Anal., 41 (2021), 1522–1545. https://doi.org/10.1093/imanum/draa013 doi: 10.1093/imanum/draa013
    [5] Y. Zhao, C. Shen, M. Qu, W. P. Bu, Y. F. Tang, Finite element methods for fractional diffusion equations, Int. J. Model., Simul., Sci. Comput., 11 (2020), 2030001. https://doi.org/10.1142/S1793962320300010 doi: 10.1142/S1793962320300010
    [6] X. Li, H. Rui, Two temporal second-order H1-Galerkin mixed finite element schemes for distributed-order fractional sub-diffusion equations, Numerical Algorithms, 79 (2018), 1107–1130. https://doi.org/10.1007/s11075-018-0476-4 doi: 10.1007/s11075-018-0476-4
    [7] W. Bu, A. Xiao, W. Zeng, Finite difference/finite element methods for distributed-order time fractional diffusion equations, J. Sci. Comput., 72 (2017), 422–441. https://doi.org/10.1007/s10915-017-0360-8 doi: 10.1007/s10915-017-0360-8
    [8] L. Feng, P. Zhuang, F. Liu, I. Turner, Y. Gu, Finite element method for space-time fractional diffusion equation, Numerical Algorithms, 72 (2016), 749–767. https://doi.org/10.1007/s11075-015-0065-8 doi: 10.1007/s11075-015-0065-8
    [9] A. J. Cheng, H. Wang, K. X. Wang, A Eulerian-Lagrangian control volume method for solute transport with anomalous diffusion, Numer. Methods Partial Differ. Equ., 31 (2015), 253–267. https://doi.org/10.1002/num.21901 doi: 10.1002/num.21901
    [10] M. Badr, A. Yazdani, H. Jafari, Stability of a finite volume element method for the time-fractional advection-diffusion equation, Numer. Methods Partial Differ. Equ., 34 (2018), 1459–1471. https://doi.org/10.1002/num.22243 doi: 10.1002/num.22243
    [11] F. Liu, P. Zhuang, I. Turner, K. Burrage, V. Anh, A new fractional finite volume method for solving the fractional diffusion equation, Appl. Math. Modell., 38 (2014), 3871–3878. https://doi.org/10.1016/j.apm.2013.10.007 doi: 10.1016/j.apm.2013.10.007
    [12] M. Stynes, E. O'Riordan, J. Gracia, Error analysis of a finite difference method on graded meshes for a time-fractional diffusion equation, SIAM J. Numer. Anal., 55 (2017), 1057–1079. https://doi.org/10.1137/16M1082329 doi: 10.1137/16M1082329
    [13] R. Choudhary, S. Singh, P. Das, D. Kumar, A higher-order stable numerical approximation for time-fractional non-linear Kuramoto-Sivashinsky equation based on quintic B-spline, Math. Methods Appl. Sci., 47 (2024), 1–23. https://doi.org/10.1002/mma.9778 doi: 10.1002/mma.9778
    [14] X. M. Gu, H. W. Sun, Y. L. Zhao, X. C. Zheng, An implicit difference scheme for time-fractional diffusion equations with a time-invariant type variable order, Appl. Math. Lett., 120 (2021), 107270. https://doi.org/10.1016/j.aml.2021.107270 doi: 10.1016/j.aml.2021.107270
    [15] X. D. Zhang, Y. L. Feng, Z. Y. Luo, J. Liu, A spatial sixth-order numerical scheme for solving fractional partial differential equation, Appl. Math. Lett., 159 (2025), 109265. https://doi.org/10.1016/j.aml.2024.109265 doi: 10.1016/j.aml.2024.109265
    [16] Y. Feng, X. Zhang, Y. Chen, L. Wei, A compact finite difference scheme for solving fractional Black-Scholes option pricing model, J. Inequal. Appl., 36 (2025), 36. https://doi.org/10.1186/s13660-025-03261-2 doi: 10.1186/s13660-025-03261-2
    [17] C. Li, F. Zeng, F. Liu, Spectral approximations to the fractional integral and derivative, Fract. Calc. Appl. Anal., 15 (2012), 383–406. doi.org/10.2478/s13540-012-0028-x doi: 10.2478/s13540-012-0028-x
    [18] S. Guo, L. Mei, Z. Zhang, Y. Jiang, Finite difference/spectral-Galerkin method for a two-dimensional distributed-order time-space fractional reaction-diffusion equation, Appl. Math. Lett., 85 (2018), 157–163. https://doi.org/10.1016/j.aml.2018.06.005 doi: 10.1016/j.aml.2018.06.005
    [19] X. Li, C. Xu, A space-time spectral method for the time fractional diffusion equation, SIAM J. Numer. Anal., 47 (2009), 2108–2131. https://doi.org/10.1137/080718942 doi: 10.1137/080718942
    [20] A. Bhardwaj, A. Kumar, A meshless method for time fractional nonlinear mixed diffusion and diffusion-wave equation, Appl. Numer. Math., 160 (2021), 146–165. https://doi.org/10.1016/j.apnum.2020.09.019 doi: 10.1016/j.apnum.2020.09.019
    [21] Y. Gu, H. G. Sun, A meshless method for solving three-dimensional time fractional diffusion equation with variable-order derivatives, Appl. Math. Modell., 78 (2020), 539–549. https://doi.org/10.1016/j.apm.2019.09.055 doi: 10.1016/j.apm.2019.09.055
    [22] V. R. Hosseini, E. Shivanian, W. Chen, Local integration of 2-D fractional telegraph equation via local radial point interpolant approximation, Eur. Phys. J. Plus, 130 (2015), 33. https://doi.org/10.1140/epjp/i2015-15033-5 doi: 10.1140/epjp/i2015-15033-5
    [23] Z. Avazzadeh, W. Chen, V. R. Hosseini, The coupling of RBF and FDM for solving higher order fractional partial differential equations, Appl. Mech. Mater., 598 (2014), 409–413. https://doi.org/10.4028/www.scientific.net/AMM.598.409 doi: 10.4028/www.scientific.net/AMM.598.409
    [24] P. Das, S. Rana, H. Ramos, A perturbation-based approach for solving fractional-order Volterra-Fredholm integro-differential equations and its convergence analysis, Int. J. Comput. Math., 97 (2020), 1994–2014. https://doi.org/10.1080/00207160.2019.1673892 doi: 10.1080/00207160.2019.1673892
    [25] P. Das, S. Rana, H. Ramos, Homotopy perturbation method for solving Caputo-type fractional-order Volterra-Fredholm integro-differential equations, Comput. Math. Methods, 1 (2019), e1047. https://doi.org/10.1002/cmm4.1047 doi: 10.1002/cmm4.1047
    [26] P. Das, S. Rana, H. Ramos, On the approximate solutions of a class of fractional order nonlinear Volterra integro-differential initial value problems and boundary value problems of first kind and their convergence analysis, J. Comput. Appl. Math., 404 (2022), 113116. https://doi.org/10.1016/j.cam.2020.113116 doi: 10.1016/j.cam.2020.113116
    [27] L. Wei, Y. F. Yang, Optimal order finite difference/local discontinuous Galerkin method for variable-order time-fractional diffusion equation, J. Comput. Appl. Math., 383 (2021), 113129. https://doi.org/10.1016/j.cam.2020.113129 doi: 10.1016/j.cam.2020.113129
    [28] L. Wei, W. Li, Local discontinuous Galerkin approximations to variable-order time-fractional diffusion model based on the Caputo-Fabrizio fractional derivative, Math. Comput. Simul., 188 (2021), 280–290. https://doi.org/10.1016/j.matcom.2021.04.001 doi: 10.1016/j.matcom.2021.04.001
    [29] W. Li, L. Wei, Analysis of Local Discontinuous Galerkin Method for the Variable-order Subdiffusion Equation with the Caputo-Hadamard Derivative, Taiwanese J. Math., 28 (2024), 1095–1110. https://doi.org/10.11650/tjm/240801 doi: 10.11650/tjm/240801
    [30] Y. Liu, M. Zhang, H. Li, J. Li, High-order local discontinuous Galerkin method combined with WSGD-approximation for a fractional sub-diffusion equation, Comput. Math. Appl., 73 (2017), 1298–1314. https://doi.org/10.1016/j.camwa.2016.08.015 doi: 10.1016/j.camwa.2016.08.015
    [31] B. P. Moghaddam, J. A. T. Machado, Extended algorithms for approximating variable order fractional derivatives with applications, J. Sci. Comput., 71 (2016), 1351–1374. https://doi.org/10.1007/s10915-016-0343-1 doi: 10.1007/s10915-016-0343-1
    [32] L. Ramirez, C. Coimbra, On the selection and meaning of variable order operators for dynamic modeling, Int. J. Differ. Equ., 2010 (2010), 1–16. https://doi.org/10.1155/2010/846107 doi: 10.1155/2010/846107
    [33] Z. Chen, J. Z. Qian, H. B. Zhan, L. W. Chen, S. H. Luo, Mobile-immobile model of solute transport through porous and fractured media, IAHS Publ., 341 (2011), 154–158.
    [34] R. Schumer, D. A. Benson, M. M. Meerschaert, B. Baeumer, Fractal mobile/immobile solute transport, Water Resour. Res., 39 (2003), 1–12. https://doi.org/10.1029/2003WR002141 doi: 10.1029/2003WR002141
    [35] Y. Zhang, D. A. Benson, D. M. Reeves, Time and space nonlocalities underlying fractional-derivative models: Distinction and literature review of field applications, Adv. Water Resour., 32 (2009), 561–581.
    [36] K. Sadri, H. Aminikhah, An efficient numerical method for solving a class of variable-order fractional mobile-immobile advection-dispersion equations and its convergence analysis, Chaos, Solitons Fractals, 146 (2021), 110896. https://doi.org/10.1016/j.chaos.2021.110896 doi: 10.1016/j.chaos.2021.110896
    [37] H. Ma, Y. Yang, Jacobi spectral collocation method for the time variable-order fractional mobile-immobile advection-dispersion solute transport model, East Asian J. Appl. Math., 6 (2016), 337–352. https://doi.org/10.4208/eajam.141115.060616a doi: 10.4208/eajam.141115.060616a
    [38] Z. G. Liu, A. J. Cheng, X. L. Li, A second order finite dfference scheme for quasilinear time fractional parabolic equation based on new fractional derivative, Int. J. Comput. Math., 95 (2017), 396–411. https://doi.org/10.1080/00207160.2017.1290434 doi: 10.1080/00207160.2017.1290434
    [39] A. Golbabai, O. Nikan, T. Nikazad, Numerical investigation of the time fractional mobile-immobile advection-dispersion model arising from solute transport in porous media, Int. J. Appl. Comput. Math., 5 (2019), 1–22. https://doi.org/10.1007/s40819-019-0635-x doi: 10.1007/s40819-019-0635-x
    [40] H. Zhang, F. Liu, M. S. Phanikumar, M. M. Meerschaert, A novel numerical method for the time variable fractional order mobile-immobile advection-dispersion model, Comput. Math. Appl., 66 (2013), 693–701. https://doi.org/10.1016/j.camwa.2013.01.031 doi: 10.1016/j.camwa.2013.01.031
    [41] W. Jiang, N. Liu, A numerical method for solving the time variable fractional order mobile-immobile advection-dispersion model, Appl. Numer. Math., 119 (2017), 18–32. https://doi.org/10.1016/j.apnum.2017.03.014 doi: 10.1016/j.apnum.2017.03.014
    [42] M. Saffarian, A. Mohebbi, An efficient numerical method for the solution of 2D variable order time fractional mobile-immobile advection-dispersion model, Math. Methods Appl. Sci., 44 (2021), 5908–5929. https://doi.org/10.1002/mma.7158 doi: 10.1002/mma.7158
    [43] C. F. M. Coimbra, Mechanica with variable-order differential operators, Ann. Phys., 12 (2003), 692–703. https://doi.org/10.1002/andp.200351511-1203 doi: 10.1002/andp.200351511-1203
    [44] Y. Cheng, X. Meng, Q. Zhang, Application of generalized Gauss-Radau projections for the local discontinuous Galerkin method for linear convection-diffusion equations, Math. Comp., 86 (2017), 1233–1267. https://doi.org/10.1090/mcom/3141 doi: 10.1090/mcom/3141
    [45] L. Wei, H. Wang, Y. Chen, Local discontinuous Galerkin method for a hidden-memory variable order reaction-diffusion equation, J. Appl. Math. Comput., 69 (2023), 2857–2872. https://doi.org/10.1007/s12190-023-01865-9 doi: 10.1007/s12190-023-01865-9
    [46] W. Wang, E. Barkai, Fractional advection-diffusion-asymmetry equation, Phys. Rev. Lett., 125 (2020), 240606. https://doi.org/10.1103/PhysRevLett.125.240606 doi: 10.1103/PhysRevLett.125.240606
  • Reader Comments
  • © 2025 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Metrics

Article views(312) PDF downloads(27) Cited by(0)

Figures and Tables

Figures(2)  /  Tables(3)

Other Articles By Authors

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog