Variable-Order Fractional Wave Equation

• Variable-order fractional wave equations are partial differential equations with derivatives whose orders vary in space and time, capturing evolving memory effects in heterogeneous media.
• They employ Caputo and Riemann–Liouville formulations integrated with advanced discretization methods such as spectral and finite element schemes for precise numerical simulation.
• Applications include modeling diffusive waves in viscoelastic media and developing fast, stable solvers using techniques like Toeplitz structures and Fourier matrix-free methods.

A variable-order fractional wave equation is a class of partial differential equations where the order of the fractional derivative—temporal or spatial—depends on position and/or time, and typically models processes with evolving memory or heterogeneity, such as diffusive waves in viscoelastic media with changing physical properties. These equations generalize classic wave and diffusion equations by allowing the fractional order $\alpha$, $\mu$, or $s$ to vary (e.g., $\alpha(t)$ or $s(x)$), introducing multi-scale, history-dependent dynamics not present in constant-order models. This article synthesizes recent theoretical, analytical, and numerical developments for variable-order fractional wave equations, including damping and diffusion-wave constructions, rigorous regularity theory, discretization methods, and high-performance algorithms.

1. Mathematical Formulations and Variable-Order Derivatives

Variable-order fractional wave equations encompass several distinct formulations:

• Variable-order time-fractional wave equation (damping):
  • Van Bockstal–Zaky–Hendy (Bockstal et al., 2023) analyze the PDE:

  $$\Phi_{tt}(x,t) + \rho(x,t) \partial_t^{\mu(t)} \Phi(x,t) = \nabla \cdot [\beta(x,t) \nabla \Phi(x,t)] + f(\Phi(x,t)) + Q(x,t),$$

  where $\partial_t^{\mu(t)}$ is the Caputo derivative of time-varying order $\mu(t)$.

  $$\partial_t^{\mu(t)} \Phi(t) = \frac{1}{\Gamma(1-\mu(t))} \int_0^t (t-s)^{-\mu(t)} \Phi'(s)\,ds,$$

  with $0 \leq \mu(t) < 1$.

• Variable-order space-fractional wave equation:

  • Zhang–Zhao–Zhou (Zhang et al., 2023) paper

  $$\partial_{tt}u(x,t) = -\kappa\,(-\Delta)^{s(x)}\,u(x,t) + f(u(x,t)), \quad x \in \Omega,$$

  where $(-\Delta)^{s(x)}$ is a fractional Laplacian with spatially variable order $s(x)$, defined in Fourier space as

  $$(-\Delta)^{s(x)} u(x) = \int_{\mathbb{R}^d} \widehat{u}(\boldsymbol\xi) |\boldsymbol\xi|^{2 s(x)} e^{2 \pi i x \cdot \boldsymbol\xi} d\boldsymbol\xi.$$

• Variable-exponent fractional diffusion-wave equation:

  • Li et al. (Zheng et al., 5 Jun 2024), Chen et al. (Jia et al., 8 Nov 2025):

  $^c_t^{\alpha(t)}u(x,t) - \Delta u(x,t) = f(x,t),\quad ^c_t^{\alpha(t)}u := (k* \partial_t^2 u)(t),$

  with $\alpha(t) \in (1,2)$ and convolution kernel $k(t) = t^{1-\alpha(t)}/\Gamma(2-\alpha(t))$.

• Variable-order Helmholtz equation (stationary wave):

  • Zhao–Mao–Karniadakis (Zhao et al., 2018):

  $\lambda^2 u(x) - D_x^{\alpha(x)} u(x) = f(x),$

  with variable-order Riemann–Liouville derivative $D_x^{\alpha(x)}$.

Variable-order derivatives are typically defined via Caputo or Riemann–Liouville formulations, but in space fractional cases, the pseudo-differential (Fourier multipliers) or singular integral kernel definitions arise. Orders can vary in time, space, or both; for $\alpha(t)$ or $s(x)$, order regularity is crucial for well-posedness and high-order error estimates.

2. Well-posedness, Regularity, and Analytical Reformulation

• Well-posedness: Existence and uniqueness require energy estimates accompanying regularity assumptions on the variable order.

  • For time-fractional damping with Caputo derivatives, Van Bockstal et al. (Bockstal et al., 2023) prove existence and uniqueness under mild boundedness conditions on $\mu(t) \in L^\infty$, coefficients $\rho, \beta$, and Lipschitz nonlinearity $f$.
  • For variable exponent models with convolution kernels (e.g., (Zheng et al., 5 Jun 2024, Jia et al., 8 Nov 2025)), regularity in $\alpha$ or $s$ is essential. Chen et al. (Jia et al., 8 Nov 2025) require $\alpha \in C^2[0,T]$, $0 \leq \alpha(t) < 1$, and eliminate weak singularity at $t=0$ via $\alpha(0) = \alpha'(0) = 0$.
• Analytical reformulation: To facilitate numerical analysis, convolution-based reformulations are introduced:
  • By convolving Abel kernels and shifting variable exponent kernels to more tractable forms (e.g., (Zheng et al., 5 Jun 2024, Jia et al., 8 Nov 2025)), well-posedness is reduced to energy and Volterra-type integro-differential equations with bounded convolution coefficients.
  • Modal decomposition in an eigenbasis (see (Jia et al., 8 Nov 2025)) yields scalar ODEs with variable-order convolution terms, enabling sharp regularity and stability results:

  $$u_i'' + \lambda_i^2 u_i = f_i - \lambda_i^2 (k * u_i).$$

• Regularity results: High-order regularity is established when initial data, forcing terms, and the fractional order are sufficiently smooth. Chen et al. (Jia et al., 8 Nov 2025) prove $u \in H^k(L^2) \cap H^{k-2}(\dot{H}^2)$ provided $u_0 \in \dot{H}^k$, $\tilde u_0 \in \dot{H}^{k-1}$, and $f \in H^{k-1}(\dot{H}^1) \cap L^2(\dot{H}^k)$.

3. Discretization: Time-stepping and Spatial Schemes

• Time discretization:

  • Rothe's method (Bockstal et al., 2023) uses backward Euler for time derivatives and discrete convolution quadrature (CQ) for Caputo damping. It yields time-discrete variational problems with explicit formulas for the discrete variable-order Caputo convolution, assembling weights $b^i_q$ for efficient recursion and stability.
  • High-order schemes (Zheng et al., 5 Jun 2024) employ BDF2 or second-order product-integration rules for variable-exponent equations, ensuring accuracy up to $\mathcal{O}(\tau^{\alpha_0})$ and $\mathcal{O}(\tau^2)$ in time under appropriate regularity.
  • Finite element Ritz–Volterra projection (Jia et al., 8 Nov 2025) discretizes both space and the Volterra time-convolution, producing schemes with explicit Toeplitz block structure, exploited for fast algorithms.
• Spatial discretization:
  • Galerkin spectral methods (Bockstal et al., 2023) utilize Legendre polynomials for $H^1_0$ projection, attaining spectral accuracy in space for smooth solutions. The discrete system at each time step is small, leading to efficient high-fidelity computation.
  • Multi-domain spectral collocation (Zhao et al., 2018) constructs penalized global differentiation matrices via Jacobi polynomial expansions, supporting $hp$-refinement and sharp layers for variable-order Riemann–Liouville operators.
  • Fourier pseudospectral algorithms (Zhang et al., 2023) represent the spatial variable-order fractional Laplacian via Taylor expansions of the Fourier multiplier, enabling spectral accuracy in space and fast matrix-free evaluation by $M\ll N$ inverse FFTs per time step.

4. Fast Algorithms and Computational Considerations

• Toeplitz structure and divide-and-conquer solvers:
  • Discretizations with convolution terms in time give block-Toeplitz systems (Jia et al., 8 Nov 2025), with memory requirements scaling with $O(MN^2)$, where $M$ is the spatial DOF and $N$ is the number of time steps.
  • A fast divide-and-conquer algorithm reduces solver complexity to $O(MN\log^2 N)$ by recursively partitioning the time domain, utilizing FFTs for block matrix multiplication in the convolution summation (see Section 8 in (Jia et al., 8 Nov 2025)).
• Accelerated matrix-free spatial operators:
  • For spatially variable fractional orders, constructing and storing the full operator matrix is computationally prohibitive ($O(J^2)$ storage for $J$ points).
  • The matrix-free approach (Zhang et al., 2023) applies local Taylor expansions in the order $s(x)$, requiring only $O(MN\log N)$ time and $O(MN)$ storage.
  • Parallelization is supported via distribution of FFTs over expansion terms.
• Penalty stabilization:
  • In multi-domain spectral collocation (Zhao et al., 2018), penalty parameters enforce flux continuity at element interfaces, regularizing the eigenvalues of the differentiation matrix and improving stability in time-dependent problems.
  • Eigenvalue monitoring and condition number control are crucial for setting penalty levels to maintain robustness without ill-conditioning.

Discretization	Key techniques	Complexity per step
Rothe–Galerkin	CQ, Legendre spectral	Small (per time step)
Ritz–Volterra FE	Toeplitz, backward Euler	$O(MN\log^2 N)$ (fast solver)
Fourier Matrix-Free	Taylor–FFT, expansion	$O(MN\log N)$ (space)
MD-Spectral Collocation	Jacobi, penalty, $hp$	Problem-dependent; $O(N^p)$

5. Error Estimates, Stability, and Convergence Results

• Error bounds and rates:
  • For Rothe–Galerkin (Bockstal et al., 2023), the spatial error is spectral ($C N^{-m}$), temporal convergence is $O(\tau^{1-\max \mu})$; degradation in space when time error dominates.
  • In high-order time discretizations (Zheng et al., 5 Jun 2024, Jia et al., 8 Nov 2025), error bounds of $O(\tau^{\alpha_0} + h^2)$ (first order in time, second order in space) are proved; second-order time-stepping achieves $O(\tau^2 + h^2)$.
  • Matrix-free Fourier methods (Zhang et al., 2023) demonstrate spectral accuracy in space and second-order accuracy in time, independent of variable order, via rigorous numerical experimentation.
• Stability:
  • Discrete energy methods (Grönwall, Young’s inequalities) are systematically invoked to guarantee stability, as shown in Theorems 3.2, 3.6 (Zheng et al., 5 Jun 2024, Jia et al., 8 Nov 2025).
  • Penalty methods in spatial schemes (Zhao et al., 2018) enforce the necessary stability conditions—failure to penalize results in positive eigenvalues and instability.
• Convergence observations:
  • For smooth solutions and sufficiently regular variable order, all schemes exhibit the designed rates (spectral or algebraic, per method).
  • Numerical evidence (Zheng et al., 5 Jun 2024) suggests that "zero-derivative" constraints on the variable order at $t=0$ may be purely technical, as convergence remains robust when these are violated.

6. Applications, Modeling Implications, and Outlook

• Modeling heterogeneous media: Variable-order fractional wave equations directly model diffusive or dispersive waves in viscoelastic media where local physical properties (memory, elasticity, attenuation) evolve in space or time.
  • Comparative studies (Zhang et al., 2023) show wavefronts in fractional heterogeneous media exhibit more intricate dynamics than homogeneous counterparts—nonlocality and local order produce complex propagation, deformation, and soliton interaction profiles.
• Practical test problems:
  • Manufactured solutions (e.g., $t^2 \sin(\pi x)$, $t^3 \sin(2\pi x)$, spatial polynomial profiles) validate numerical convergence and stability.
  • Oscillatory, linear, and quadratic forms for variable order ($\mu(t)$ or $\alpha(t)$) are empirically studied (Bockstal et al., 2023).
• Numerical efficiency: Fast algorithms and matrix-free strategies render high-dimensional computations feasible, as direct approaches are prohibitive for large $N$ or variable $s(x)$.
• Remaining questions and opportunities:
  • Extension to fully variable order in both space and time, nonhomogeneous boundary conditions, and generalized nonlocal operators are designated as open research areas (Zheng et al., 5 Jun 2024).
  • Removal of technical regularity constraints (e.g., at $t=0$) and deeper error theory for singular solutions are avenues for development.
• Impact: These advances establish a robust, high-order accurate methodology for variable-order fractional wave modeling, bridging rigorous theory and scalable computational practice for heterogeneous, evolving media in applied mathematics, geophysics, and material science.