Geometry-Aware Stability and Stiffness-Robust Fractional Time Integration of Partial Differential Equations by Dual-Reciprocity Boundary Elements and Implicit Adams–Moulton Schemes

1. Introduction

Fractional partial differential equations have become a mature modeling paradigm for anomalous diffusion, hereditary transport, relaxation, thermal memory, viscoelastic response, porous-media dynamics, and multiscale processes in which the present state depends on a nonlocal temporal history. Their numerical solution is substantially more demanding than that of classical parabolic equations because every new time level is coupled to the entire preceding trajectory and because the exact solution commonly possesses limited regularity near the initial state. The space-spectral time Adams–Bashforth–Moulton formulation of Sohail, Maqbool, and Ellahi provided an important stability-oriented reference for fractional partial differential equations [1]. The modern fractional Adams methodology is supported by detailed predictor–corrector error analysis [2], stability results for fractional Adams–Moulton formulas [3], and explicit characterization of predictor–corrector stability regions [4]. Graded-mesh analysis has clarified the order degradation induced by weak initial singularities [5], while fast high-order predictor–corrector constructions have demonstrated that graded temporal resolution and compressed memory can be combined without abandoning the Adams framework [6]. These developments establish a strong temporal foundation, but they do not by themselves resolve the interaction between fractional time integration and a stiff spatial operator generated on a geometrically complex domain.

For irregular, multiply connected, curved, or nonconvex domains, global spectral collocation may become restrictive because geometric mappings, global basis construction, and boundary-condition enforcement can dominate the numerical design. Boundary element methods offer a fundamentally different spatial architecture by concentrating geometric discretization on the boundary while retaining direct access to boundary fluxes and interior reconstruction. The dual-reciprocity formulation of Wrobel and Brebbia transformed nonhomogeneous diffusion terms into boundary quantities [7], and early nonlinear transient applications showed that the same mechanism could accommodate strongly coupled thermal source terms [8]. More recent fractional boundary element developments have treated Caputo anomalous diffusion directly [9], extended boundary integral formulations to fractional diffusion-wave dynamics [10], and combined boundary reduction with global temporal algebra for multiorder advection–diffusion–reaction equations [11]. Collectively, these studies demonstrate the geometric flexibility of boundary formulations; however, they also expose a decisive methodological question that is often left implicit: the time integrator must be designed for the actual dense and potentially nonnormal evolution operator produced by the boundary discretization.

A second difficulty is the hereditary memory cost. Direct evaluation of the Caputo history produces growing storage and arithmetic burdens as the simulation horizon increases. Sum-of-exponentials acceleration provides a controlled mechanism for reducing this cost while preserving fractional accuracy [12]. At the same time, explicit fractional diffusion schemes can possess stability restrictions that are substantially more severe than their integer-order analogues [13]. Convolution quadrature provides an alternative operator-based route to stable fractional evolution [14], founded on Lubich’s discretized operational calculus [15], while adaptive multidimensional simulations have emphasized the need to balance temporal history resolution against spatial approximation error [16]. The resulting conclusion is structural rather than algorithm-specific: memory, stiffness, and spatial discretization are coupled features of the same fully discrete problem and should therefore be analyzed within one operator-level framework.

Temporal regularity further complicates the design of high-order algorithms. The Alikhanov difference formula provides a widely used stable high-order benchmark for time-fractional diffusion [17]. Sharp graded-mesh analysis quantifies the interaction between start-up singularity and temporal refinement [18], while nonuniform first-order theory [19] and discrete fractional Grönwall inequalities [20] provide powerful stability tools for variable-step discretizations. Earlier fully discrete diffusion-wave analysis [21], finite-difference spectral approximations [22], higher-order Caputo differentiation formulas [23], and time-stepping estimates for nonsmooth initial data [24] all show that formal consistency order is not sufficient to guarantee practical reliability. A credible fractional solver must state its regularity assumptions explicitly, separate smooth compatible verification from generic weakly singular behavior, and provide a mechanism for recovering accuracy when the initial layer is unresolved. The overview of Jin, Lazarov, and Zhou situates these issues within the broader numerical theory of time-fractional evolution equations with nonsmooth data [25].

Transient boundary element research has addressed dense algebra and time marching through several complementary strategies. Matrix inversion can be avoided or reorganized in dual-reciprocity heat conduction [26], precise integration can be applied to the semidiscrete dual-reciprocity system [27], and isogeometric dual reciprocity can preserve exact geometric descriptions while handling transient heat sources [28]. Earlier nonlinear boundary element formulations already demonstrated the effectiveness of time-independent fundamental solutions for transient nonlinear heat conduction [29]. More recently, nonuniform second-order fractional formulas have supplied an additional route for resolving early-time singular behavior without sacrificing asymptotic accuracy [30]. These results strongly motivate a fractional boundary element formulation in which geometry reduction, stiffness-aware time integration, memory acceleration, and post-assembly stability certification are treated as mutually dependent components rather than as sequential implementation choices.

Fractional multistep analysis provides the appropriate temporal machinery for such a construction. General fractional multistep formulas and their stability properties were developed by Galeone and Garrappa [31], whereas the Diethelm–Ford–Freed predictor–corrector method remains a foundational practical scheme [32]. A broader collection of fractional numerical algorithms was subsequently assembled by Diethelm and coauthors [33], and Garrappa’s survey clarified implementation choices among product-integration and fractional linear multistep families [34]. Best-practice guidance has further emphasized that intuition inherited from integer-order time stepping can be misleading for fractional operators [35]. This caution is especially important for dense boundary element evolution matrices because collocation, mixed boundary conditions, advection, and geometric asymmetry can produce nonnormal operators whose transient behavior is not characterized adequately by eigenvalues alone. Pseudospectral analysis supplies the natural language for such transient amplification [36], while numerical-range estimates provide a complementary functional-calculus control mechanism [37].

The dual-reciprocity interpolation stage introduces its own approximation and conditioning mechanisms. Multiquadric radial basis constructions established highly accurate scattered-data interpolation and derivative recovery [38] and demonstrated their effectiveness for parabolic, hyperbolic, and elliptic partial differential equations [39]. Compactly supported positive-definite radial functions subsequently provided a systematic mechanism for improving locality and conditioning [40]. These developments support radial basis interpolation inside the dual-reciprocity transformation while also indicating clear routes toward scalable variants based on localized interpolation, compact support, adaptive point placement, and matrix compression.

In this paper, a stability-centered dual-reciprocity boundary element Adams–Moulton framework is developed for time-fractional diffusion–reaction equations. The central departure from the space-spectral Adams formulation is not merely the replacement of one spatial discretization by another. Instead, the stiff boundary element diffusion operator is retained implicitly in the Moulton correction, whereas nonlinear or lower-stiffness contributions are evaluated through an Adams predictor. The resulting semidiscrete fractional system is analyzed through the spectrum, numerical range, and resolvent of the assembled boundary element evolution operator. This geometry-aware viewpoint converts stability from an abstract scalar property into a computable post-assembly certificate that reflects the actual domain, boundary conditions, collocation structure, and operator nonnormality. It also removes the dominant explicit diffusive restriction from the predictor–corrector architecture while preserving a boundary-reduced geometric representation.

The contributions of the paper are therefore organized around a unified analytical and computational objective. First, a complete dual-reciprocity reduction is derived for a Caputo diffusion–reaction model and expressed in canonical matrix fractional evolution form. Second, a linearly implicit Adams correction is constructed so that the dense stiff operator is factorized once for fixed geometry and time increment, while nonlinear terms remain compatible with predictor evaluation and history compression. Third, modal, numerical-range, and resolvent-based stability statements are developed together with a combined space–time convergence description that distinguishes smooth compatible solutions from weakly singular initial layers. Fourth, a practical post-assembly stability workflow is formulated to connect the theoretical stability region directly to the computed boundary element operator. Fifth, the method is verified on a curved disk and a nonconvex L-shaped domain through manufactured-solution accuracy, spatial and temporal convergence, quantified stiffness growth, spectrum-based diagnostics, memory-retention behavior, and comparison with explicit fractional predictor–corrector limits. The resulting framework advances the earlier stability analysis from a geometry-independent setting to an operator-level boundary element methodology suitable for complex domains and extensible multiphysics models.

2. Model Class and Fractional Preliminaries

Let $\Omega \subset {\mathbb{R}}^2$ be a bounded Lipschitz domain with its boundary $\Gamma ={\Gamma }_D\cup {\Gamma }_N$. Then

$$D_t^{\alpha ,\mathrm{C}}u\left(\mathbf{x},t\right)=\kappa \Delta u\left(\mathbf{x},t\right)+\mathcal{R}\left(\mathbf{x},t,u,\nabla u\right)+s\left(\mathbf{x},t\right), \mathbf{x}\in \Omega , 0<t\le T,$$

(1)

where $0<\alpha <1$, $\kappa >0$, and $\mathcal{R}$ may include reaction, advection, weak nonlinearities, or additional constitutive source terms. The initial and mixed boundary data are

$$u\left(\mathbf{x},0\right)=u_0\left(\mathbf{x}\right), u=g_D\hspace{0.25em}\mathrm{on} {\Gamma }_D, q={\displaystyle \frac{\partial u}{\partial n}}=g_N\hspace{0.25em}\mathrm{on} {\Gamma }_N.$$

(2)

For $u$ sufficiently regular, the Caputo derivative is

$$D_t^{\alpha ,\mathrm{C}}u\left(t\right)={\displaystyle \frac{1}{\Gamma \left(1-\alpha \right)}}{\int }_0^t{\left(t-s\right)}^{-\alpha }u\mathrm{\text{'}}\left(s\right)\hspace{0.17em}ds.$$

(3)

Equation (1) is equivalent to the Volterra integral equation

$$u\left(t\right)=u_0+{\displaystyle \frac{1}{\Gamma \left(\alpha \right)}}{\int }_0^t{\left(t-s\right)}^{\alpha -1}\left[\kappa \Delta u\left(s\right)+\mathcal{R}\left(s,u\left(s\right)\right)+s\left(s\right)\right]ds.$$

(4)

The fractional Adams family approximates the history integral in (4) by product integration. The difficulty is that the right-hand side after spatial discretization contains a stiff matrix. In the present construction that matrix is generated by DRBEM and is treated implicitly in the corrector.

3. Dual-Reciprocity Boundary Element Spatial Reduction

3.1. Boundary Integral Identity

Rewrite (1) as

$$\Delta u=b\left(\mathbf{x},t\right), b={\kappa }^{-1}\left(D_t^{\alpha ,\mathrm{C}}u-\mathcal{R}-s\right).$$

(5)

For the two-dimensional Laplacian, take the fundamental solution

$$u^{\mathrm{*}}\left(\mathbf{x},\mathsf{\xi }\right)=-{\displaystyle \frac{1}{2\pi }}\ln r, r=\parallel \mathbf{x}-\mathsf{\xi }\parallel .$$

(6)

Green’s identity gives, for a source point $\mathsf{\xi }$,

$$c\left(\mathsf{\xi }\right)u\left(\mathsf{\xi }\right)+{\int }_{\Gamma }u\hspace{0.17em}{q}^{\mathrm{*}}\hspace{0.17em}d\Gamma -{\int }_{\Gamma }q\hspace{0.17em}{u}^{\mathrm{*}}\hspace{0.17em}d\Gamma =\mathcal{D}\left[b\right]\left(\mathsf{\xi }\right),$$

(7)

where $c=1/2$ at a smooth boundary point and $c=1$ in the interior, and $\mathcal{D}\left[b\right]$ denotes the domain contribution with the sign convention induced by (6).

A direct domain BEM would require cells to evaluate $\mathcal{D}\left[b\right]$. Dual reciprocity avoids this by approximating $b$ at boundary and selected interior interpolation points.

3.2. Dual Reciprocity Approximation

Let $\{{\mathbf{x}}_j{\}}_{j=1}^{N_c}$ contain all boundary collocation points and $N_i$ interior dual-reciprocity points. Approximate

$$b\left(\mathbf{x},t\right)\approx {\displaystyle \sum _{j=1}^{N_c}}{\beta }_j\left(t\right){\varphi }_j\left(\mathbf{x}\right), \mathsf{\beta }=F^{-1}\mathbf{b},$$

(8)

with $F_{ij}={\varphi }_j\left({\mathbf{x}}_i\right)$. The numerical experiments use the classical radial basis

$${\varphi }_j\left(\mathbf{x}\right)=1+r_j, r_j=\parallel \mathbf{x}-{\mathbf{x}}_j\parallel .$$

(9)

In two dimensions a convenient particular solution is

$${\widehat{u}}_j={\displaystyle \frac{r_j^2}{4}}+{\displaystyle \frac{r_j^3}{9}}, \Delta {\widehat{u}}_j=1+r_j={\varphi }_j.$$

(10)

Thus the domain term can be converted to boundary quantities of ${\widehat{u}}_j$ and ${\widehat{q}}_j=\partial {\widehat{u}}_j/\partial n$. After boundary discretization,

$$H{\mathbf{u}}_b-G{\mathbf{q}}_b=D\mathbf{b}, D=\left(H\widehat{U}-G\widehat{Q}\right)F^{-1},$$

(11)

with the usual interior $c=1$ contribution included in the corresponding rows of $H\widehat{U}$.

For straight constant boundary elements, all nonsingular integrals are evaluated by Gaussian quadrature. The singular diagonal of $G$ has the analytical value

$$G_{ii}={\displaystyle \frac{L_i}{2\pi }}\left[1-\ln \left({\displaystyle \frac{L_i}{2}}\right)\right],$$

(12)

where $L_i$ is the element length. Linear or quadratic boundary elements can be substituted without changing the time integration developed below.

Figure 1. The collocation architecture illustrates the boundary-only geometric discretization supplemented by interior dual-reciprocity interpolation points. The distribution is selected to resolve the nonhomogeneous fractional source field while preserving the principal BEM advantage of avoiding a conventional volume-element mesh.

3.3. Evolution Operator After Boundary-Condition Elimination

To show the algebra transparently, first assume homogeneous Dirichlet data on the full boundary. Partition rows of (11) into boundary ($b$) and interior ($i$) equations. Since $u_b=0$ and $D_t^{\alpha ,\mathrm{C}}{u}_b=0$,

$$-G_{bb}{\mathbf{q}}_b=D_{bi}\hspace{0.17em}{D}_t^{\alpha ,\mathrm{C}}\mathbf{U}-D_b\mathbf{f},$$

(13)

$$\mathbf{U}-G_{ib}{\mathbf{q}}_b=D_{ii}\hspace{0.17em}{D}_t^{\alpha ,\mathrm{C}}\mathbf{U}-D_i\mathbf{f},$$

(14)

where $\mathbf{U}$ contains interior nodal values and $\mathbf{f}$ collects the discretized $\mathcal{R}+s$ terms with the factor ${\kappa }^{-1}$ absorbed into the matrices for clarity. Eliminating ${\mathbf{q}}_b$ gives

$$T\hspace{0.17em}{D}_t^{\alpha ,\mathrm{C}}\mathbf{U}=\mathbf{U}-R\mathbf{f},$$

(15)

where

$$T=D_{ii}-G_{ib}{G}_{bb}^{-1}{D}_{bi}, R=G_{ib}{G}_{bb}^{-1}{D}_b-D_i.$$

(16)

Hence the semidiscrete equation can be written in canonical form

$$D_t^{\alpha ,\mathrm{C}}\mathbf{U}=A_h\mathbf{U}+{\mathbf{N}}_h\left(t,\mathbf{U}\right)+{\mathbf{g}}_h\left(t\right), A_h=\kappa T^{-1}.$$

(17)

For general mixed boundary conditions, prescribed and unknown boundary variables are exchanged in the standard BEM manner before the Schur complement is formed. The essential point is unchanged: the FPDE is reduced to a finite-dimensional Caputo system whose stiff linear part is the assembled boundary-element operator $A_h$.

4. DRBEM-IABM Time Integration

4.1. Fractional Adams-Bashforth Predictor

Let $t_n=n\tau$ and define

$${\mathcal{F}}_n=A_h{\mathbf{U}}_n+{\mathbf{N}}_h\left(t_n,{\mathbf{U}}_n\right)+{\mathbf{g}}_h\left(t_n\right).$$

(18)

The one-step fractional Adams-Bashforth predictor is

$${\mathbf{U}}_{n+1}^P={\mathbf{U}}_0+{\displaystyle \frac{{\tau }^{\alpha }}{\Gamma \left(\alpha +1\right)}}{\displaystyle \sum _{j=0}^n}{b}_{j,n+1}{\mathcal{F}}_j,$$

(19)

$$b_{j,n+1}={\left(n+1-j\right)}^{\alpha }-{\left(n-j\right)}^{\alpha }.$$

(20)

If (19) is followed by a corrector that evaluates the entire right-hand side only at ${\mathbf{U}}_{n+1}^P$, the resulting PECE method remains explicit with respect to $A_h$ and therefore inherits a bounded stability interval.

4.2. Linearly Implicit Fractional Adams-Moulton Corrector

The proposed modification retains the predictor for the nonlinear term but treats $A_h\mathbf{U}$ implicitly. Define

$$\omega ={\displaystyle \frac{{\tau }^{\alpha }}{\Gamma \left(\alpha +2\right)}}.$$

(21)

The corrector is

$$\left(I-\omega A_h\right){\mathbf{U}}_{n+1}={\mathbf{U}}_0+\omega \left[{\mathbf{N}}_h\left(t_{n+1},{\mathbf{U}}_{n+1}^P\right)+{\mathbf{g}}_h\left(t_{n+1}\right)+{\displaystyle \sum _{j=0}^n}{a}_{j,n+1}{\mathcal{F}}_j\right],$$

(22)

where

$$a_{0,n+1}=n^{\alpha +1}-\left(n-\alpha \right){\left(n+1\right)}^{\alpha },$$

(23a)

$$a_{j,n+1}={\left(n-j+2\right)}^{\alpha +1}+{\left(n-j\right)}^{\alpha +1}-2{\left(n-j+1\right)}^{\alpha +1}, 1\le j\le n.$$

(23b)

The coefficient multiplying the new-time right-hand side is one, as already incorporated in (22). Because $A_h$ is constant for a fixed geometry and linear diffusion operator, $I-\omega A_h$ is factorized once. Every later correction requires only triangular substitutions plus history accumulation. This is the crucial stiffness modification.

For a nonlinear diffusion coefficient, $A_h$ can be frozen over a time step, updated by a Picard/Newton outer iteration, or split into a constant stiff reference part plus an explicitly evaluated residual. The stability analysis below applies directly to the constant linear part and perturbatively to the nonlinear remainder.

4.3. Optional Graded Start-Up and Fast History Compression

Generic Caputo solutions often contain an initial layer, and a uniform fractional Adams mesh can lose its optimal order. A graded grid $t_n=T{\left(n/N\right)}^r$ can be used during a short start-up interval, following the error analysis in [5], before switching to a uniform step for constant-factorization efficiency.

Long-time simulations also face an $O\left(N_t^2\right)$ history cost. The kernel in the fractional integral can be approximated by a sum of exponentials,

$${\left(t-s\right)}^{\alpha -1}\approx {\displaystyle \sum _{\mathcal{l}=1}^{N_{\mathrm{e}\mathrm{x}\mathrm{p}}}}{w}_{\mathcal{l}}{e}^{-{\eta }_{\mathcal{l}}\left(t-s\right)},$$

(24)

which converts the distant history to a set of recursively updated auxiliary states. SOE constructions of this type reduce storage from $O\left(N_t\right)$ histories per degree of freedom to $O\left(N_{\mathrm{e}\mathrm{x}\mathrm{p}}\right)$ and the history work from quadratic to nearly linear-logarithmic complexity [16]. This acceleration is orthogonal to the implicit BEM solve and can therefore be added without changing the stability mechanism in Section 5.1, Section 5.2 and Section 5.3. The numerical validation below uses direct history accumulation so that no approximation from SOE compression is mixed into the reported convergence data.

4.4. Algorithm

Algorithm 1: DRBEM-IABM for a Caputo diffusion-reaction problem

Discretize $\Gamma$ by boundary elements and select interior DR points. Assemble $H$, $G$, $F$, $\widehat{U}$, $\widehat{Q}$ and $D=\left(H\widehat{U}-G\widehat{Q}\right)F^{-1}$. Apply boundary conditions, eliminate boundary flux/value unknowns, and form $A_h$, ${\mathbf{N}}_h$, and ${\mathbf{g}}_h$ in (17). Select $\tau$, compute $\omega$, and factorize $I-\omega A_h$ once. For $n=0,1,\dots ,N_t-1$: evaluate the Adams-Bashforth predictor (19); evaluate explicit/nonlinear terms at the predictor; accumulate the fractional Adams-Moulton history weights (23); solve (22) for ${\mathbf{U}}_{n+1}$; store or compress the new history contribution. Recover boundary fluxes and any requested interior field values from the BEM representation.

5. Stability Analysis

5.1. Scalar Test Equation and the Explicit Stiffness Barrier

Consider the fractional Dahlquist test equation

$$D_t^{\alpha ,\mathrm{C}}y=\lambda y, y\left(0\right)=1,$$

(25)

with $z=\lambda {\tau }^{\alpha }$. For the explicit PECE Adams scheme, stability is determined by a bounded region $S_{\mathrm{PECE}}^{\alpha }$ in the $z$-plane [4]. Along the negative real axis, a numerical trace of the boundary provides the practical thresholds in Table 1. The values were generated using 500 time steps and a boundedness/decay criterion; they are reported as numerical stability diagnostics rather than as replacements for the analytical stability boundary.

Table 1 summarizes the negative-real-axis boundary of the explicit fractional PECE stability region. The tabulated thresholds provide the scalar quantity required to convert each assembled BEM eigenvalue into a practical explicit time-step restriction.

For a semidiscrete eigenvalue ${\lambda }_k<0$, a necessary explicit stiffness constraint is therefore approximately

$$\tau \le {\left({\displaystyle \frac{r_{\alpha }}{\left|{\lambda }_k\right|}}\right)}^{1/\alpha }.$$

(26)

Equation (26) is particularly severe for small $\alpha$ because $1/\alpha >1$. The stiffness mechanism is therefore more restrictive than the familiar integer-order explicit diffusion condition.

Figure 2 shows a scalar example at $\alpha =0.6$ and $z=-2$. The explicit PECE solution diverges rapidly, whereas the linearly implicit Moulton correction decays.

Figure 2. The scalar stability experiment isolates the temporal stiffness mechanism before spatial coupling. The explicit predictor-corrector trajectory leaves its bounded stability region, whereas the implicit Moulton correction preserves the physically expected dissipative response under the same fractional test configuration.

5.2. Modal Stability of the BEM Semidiscretization

For the homogeneous linear problem

$$D_t^{\alpha ,\mathrm{C}}\mathbf{U}=A_h\mathbf{U},$$

(27)

suppose first that $A_h$ is diagonalizable,

$$A_h=V\Lambda V^{-1}, \Lambda =\mathrm{diag}\left({\lambda }_1,\dots ,{\lambda }_m\right).$$

(28)

The fractional Adams-Moulton corrector acts independently on each modal component. This yields the following statement.

Theorem 1

(spectral transfer of scalar stability). Let $S_{AM}^{\alpha }$ be the scalar stability set of the fractional Adams-Moulton corrector used in (22). If

$${\tau }^{\alpha }{\lambda }_k\in S_{\mathrm{AM}}^{\alpha }, k=1,\dots ,m,$$

(29)

then the fully discrete homogeneous DRBEM-IABM solution is stable, with

$$\parallel {\mathbf{U}}_n\parallel \le \kappa \left(V\right)\underset{k}{\max }\left|{\mathcal{R}}_n\left({\tau }^{\alpha }{\lambda }_k\right)\right|\hspace{0.17em}\parallel {\mathbf{U}}_0\parallel ,$$

(30)

where${\mathcal{R}}_n$is the scalar fractional discrete resolvent and$\kappa \left(V\right)=\parallel V\parallel \parallel V^{-1}\parallel$ .

Proof sketch. Transform (27) with $\mathbf{Y}=V^{-1}\mathbf{U}$. The history recurrence becomes $m$ independent scalar Adams-Moulton recurrences. Apply scalar stability to each component and transform back. $\mathrm{▫}$Galeone and Garrappa [3] proved A-stability for the relevant low-step fractional Adams-Moulton members. Consequently, for a diffusion-dominated BEM operator whose spectrum lies in the open left half-plane, the implicit corrector does not suffer from the bounded negative-axis interval in Table 1. This is the principal stability advantage of (22).

5.3. Nonnormal Operators and a Resolvent Criterion

Advection, mixed boundary conditions, geometric asymmetry, and numerical collocation can make $A_h$ nonnormal. Eigenvalues alone may then underestimate transient amplification. A more robust sufficient condition is expressed through the numerical range $W\left(A_h\right)$ and the corrector resolvent.

Theorem 2

(resolvent stability condition). Suppose there exists $\mu >0$ such that

$$\mathrm{Re}⟨A_{h}v,v⟩\le -\mu \parallel v{\parallel }^2$$

(31)

in an energy-equivalent inner product. Then

$$\parallel {\left(I-\omega A_h\right)}^{-1}\parallel \le {\displaystyle \frac{1}{1+\omega \mu }}.$$

(32)

If the explicit nonlinear term is Lipschitz with constant $L_N$ and the predictor is bounded, perturbations in (22) satisfy a discrete fractional Volterra inequality with step factor

$$\chi ={\displaystyle \frac{\omega L_N}{1+\omega \mu }}.$$

(33)

For $\chi <1$, the nonlinear correction is contractive at the new time level; the remaining history terms are controlled by a discrete fractional Gronwall estimate.

The practical message is that the stiff diffusive scale $\parallel A_h\parallel$ occurs in the stabilizing denominator of (32), not in an explicit time-step restriction. Only the explicitly treated nonlinear/nonstiff scale enters (33).

5.4. Stability Certificate Computed from the Assembled BEM Matrix

The proposed workflow therefore includes an inexpensive stability certificate after BEM assembly:

• estimate the rightmost spectral abscissa and, if necessary, the numerical range of $A_h$;
• verify that the spectrum lies in the fractional Adams-Moulton stability sector;
• estimate $\parallel {\left(I-\omega A_h\right)}^{-1}\parallel$ or a coercivity constant $\mu$;
• enforce only the nonlinear predictor restriction implied by (33).

This differs fundamentally from a geometry-independent scalar discussion: the actual discretized geometry and boundary conditions enter the stability test through $A_h$.

6. Consistency and Convergence

Let $h_B$ denote the maximum boundary-element length, and let ${\epsilon }_{\mathrm{DR}}$ denote the error of the dual-reciprocity approximation of the nonhomogeneous term. For smooth boundaries, smooth data, and a stable collocation arrangement, constant boundary elements with the radial basis (9) exhibit second-order behavior in the interior for the benchmark used in Section 7.

For the time discretization, the classical fractional Adams analysis gives the general uniform-grid order $p_t=1+\alpha$ for $0<\alpha <1$ under the minimal regularity setting [2,5]. If the compatibility conditions remove the leading fractional initial singularity and the fractional right-hand side is sufficiently smooth, the product-integration corrector exhibits the full second-order behavior seen in the present manufactured solution.

Theorem 3

(combined error form). Assume: (i) the DRBEM spatial operator is stable and second-order accurate on the selected smooth test problem; (ii) the dual-reciprocity interpolation error is bounded by ${\epsilon }_{DR}$ ; (iii) the exact semidiscrete solution satisfies the smoothness required by the fractional Adams error theorem; and (iv) the stability conditions ofSection 5 hold. Then

$$\underset{0\le n\le N_t}{\max }\parallel u\left(t_n\right)-U_n\parallel \le C\left(h_B^2+{\epsilon }_{\mathrm{DR}}+{\tau }^{p_t}\right),$$

(34)

with $p_t=1+\alpha$ generally and $p_t=2$ for the compatible smooth-data case used in Section 7. If a weak initial layer is present, an appropriately graded start-up mesh recovers the optimal rate [5].

The proof follows the standard decomposition into spatial consistency error, dual-reciprocity interpolation error, and the stable fractional Volterra time discretization error. The implicit treatment of$A_h$changes the stability constant but not the local product-integration consistency order.

7. Numerical Experiments

All results in this section were generated specifically for the present study with a two-dimensional constant-element DRBEM implementation. Sixteen-point Gauss-Legendre quadrature was used for nonsingular boundary integrals; the singular diagonal of $G$ was integrated analytically using (12). The dual-reciprocity basis was $1+r$, with particular solution (10). Dense linear algebra was used deliberately so that acceleration techniques did not obscure the stability behavior. Unless stated otherwise, $\kappa =1$.

7.1. Experiment 1: Manufactured Caputo Diffusion on the Unit Disk

Let $\Omega =\{\left(x,y\right):x^2+y^2<1\}$, with homogeneous Dirichlet boundary conditions. Choose

$$u_{\mathrm{ex}}\left(x,y,t\right)=\left(1-x^2-y^2\right)t^2.$$

(35)

Because

$$D_t^{\alpha ,\mathrm{C}}{t}^2={\displaystyle \frac{\Gamma \left(3\right)}{\Gamma \left(3-\alpha \right)}}{t}^{2-\alpha }, \Delta \left(1-x^2-y^2\right)=-4,$$

(36)

the forcing in (1) is

$$s\left(x,y,t\right)=\left(1-x^2-y^2\right){\displaystyle \frac{\Gamma \left(3\right)}{\Gamma \left(3-\alpha \right)}}{t}^{2-\alpha }+4t^2.$$

(37)

7.1.1. Spatial Convergence

The time step was fixed at $N_t=200$ for $T=1$ and $\alpha =0.6$. Table 2 reports the relative interior $L_2$ and maximum errors as the boundary discretization is refined. The number of DR points rises with the boundary resolution.

Table 2 quantifies spatial refinement of the curved-domain manufactured solution. Relative and maximum-norm errors are reported together with the observed order so that the boundary discretization and the dual-reciprocity interpolation can be assessed independently of temporal error.

The values obtained are almost two for the entire refinement procedure. Figure 3 shows the same data against an $N_b^{-2}$ reference slope.

Figure 3. The refinement curve translates the tabulated spatial errors into an asymptotic convergence diagnostic. Agreement with the reference slope confirms that the selected constant-element DRBEM and dual-reciprocity interpolation recover approximately second-order spatial behavior for the smooth curved-domain manufactured solution.

7.1.2. Temporal Convergence

To isolate temporal error from the spatial error, the $N_b=64$ discretization was held fixed and each time-step solution was compared at $T=1$ with a semidiscrete reference computed using $N_t=1280$. Table 3 gives the relative error for three fractional orders.

Table 3 isolates the temporal discretization error by holding the spatial operator fixed and comparing each time-step solution with a highly resolved semidiscrete reference. The results test robustness with respect to the fractional order and verify the compatible-data convergence regime used in the analysis.

Figure 4. Temporal refinement is evaluated against a highly resolved semidiscrete reference for several fractional orders. The nearly parallel convergence curves show that the implicit correction retains the expected high-order behavior for smooth compatible data while remaining insensitive to the stiff diffusive spectrum.

7.2. Experiment 2: Quantified Stiffness of the DRBEM Operator

Table 4 reports the magnitude of the most negative real eigenvalue of the disk DRBEM operator under refinement. The explicit time-step estimate is computed from (26) using $\alpha =0.6$ and $r_{0.6}=1.4296$ from Table 1.

Table 4 measures how refinement moves the most negative BEM eigenvalue deeper into the stiff dissipative spectrum. The associated explicit PECE estimate demonstrates that spatial refinement rapidly converts the nominally simple predictor-corrector strategy into a prohibitively small-step method.

For $N_b=64$, the implicit corrector attains a temporal error of only $3.38\times {10}^{-5}$ with $N_t=20$ in Table 3, whereas the explicit PECE estimate suggests more than ${10}^5$ steps merely to remain inside the negative-real-axis stability interval. The comparison is not intended as a CPU benchmark; it isolates the stability mechanism and shows why an implicit diffusion corrector is necessary for a practical BEM fractional solver.

7.3. Experiment 3: Nonconvex L-Shaped Domain

Consider the nonconvex domain with vertices

$$\left(-1,-1\right),\hspace{0.25em}\left(1,-1\right),\hspace{0.25em}\left(1,0\right),\hspace{0.25em}\left(0,0\right),\hspace{0.25em}\left(0,1\right),\hspace{0.25em}\left(-1,1\right),$$

(38)

with homogeneous Dirichlet boundary conditions, zero forcing, and initial state

$$u_0\left(x,y\right)=xy\left(1-x^2\right)\left(1-y^2\right).$$

(39)

The boundary is divided into 48 constant elements and 108 interior unknowns are generated by the DR point set. The assembled operator has

$$\mathrm{minRe}\lambda \left(A_h\right)=-1830.96, \mathrm{maxRe}\lambda \left(A_h\right)=-9.61, \max \left|\mathrm{Im}\lambda \left(A_h\right)\right|=56.66.$$

(40)

Thus the entire computed spectrum remains in the left half-plane despite the nonconvex geometry and nonnormal collocation effects. Figure 5 plots the spectrum.

Figure 5. The eigenvalue cloud provides a geometry-specific stability certificate for the nonconvex domain. The spectrum remains confined to the dissipative half-plane, while its spread and departure from normality motivate the complementary resolvent and numerical-range criteria developed in the stability analysis.

At $\alpha =0.6$, $T=0.5$, temporal refinement against an $N_t=800$ reference gives Table 5.

Table 5 evaluates temporal convergence on the nonconvex domain, where the re-entrant corner and nonnormal spatial operator make the test substantially more demanding than the smooth disk benchmark. The refined sequence demonstrates recovery toward the asymptotic temporal regime.

The pre-asymptotic convergence rate is influenced by the corner singularity and the nonnormal spatial operator, but the improved convergence rate exhibits the second-order temporal behavior. A representative solution at $T=0.5$ is shown in Figure 6.

Figure 6. Figure 6. The reconstructed domain shows that the boundary reduced fractional model maintains stability on the re-entrant corner geometry without using a volume mesh tailored to the domain. The solution retains the expected corner-influenced spatial structure while the fractional history produces smooth long-memory decay.

To illustrate the effect of fractional memory, the same L-shaped calculation was continued to $T=1$ with $N_t=200$. With an initial norm of ∥U₀∥₂ = 0.7704, the solution decayed to final values of 0.02834, 0.01741, and 0.00436 for α = 0.3, 0.6, and 0.9, respectively. These results show that lower fractional orders preserve a greater proportion of the initial state over time, reflecting the slower relaxation and diffusion dynamics typically associated with anomalous transport.

7.4. Cross-Benchmark Stiffness Synthesis

The refinement data reveal a systematic separation between spatial accuracy and explicit temporal feasibility. As the boundary resolution increases, the dissipative spectral radius grows and the stability-preserving explicit PECE step count increases by more than an order of magnitude over the tested sequence. This trend is not a defect of the spatial discretization; it is the expected signature of a progressively better resolved diffusion operator. The linearly implicit Moulton correction method is therefore necessary since it enables mesh refinement in space while being able to take large time steps independent of the stiffest diffusive mode.

Figure 7. The refinement curve converts the stiffness data of Table 4 into a direct computational burden metric. The logarithmic vertical scale shows the rapid growth in the number of explicit PECE steps required merely to remain inside the scalar stability interval as the boundary discretization is refined.

8. Comparison with Related Fractional BEM and Spectral Approaches

Table 6 positions the proposed method relative to representative literature.

Table 6 positions the proposed formulation against representative fractional spectral, direct BEM, hybrid BEM, and transient DRBEM strategies. The comparison emphasizes the distinctive combination of boundary-reduced geometry handling, implicit treatment of the stiff diffusion operator, and operator-level stability certification.

The contribution is therefore not simply “BEM instead of spectral collocation.” It is the coupling of boundary-reduced geometry handling with a time integrator whose stability architecture is matched to the dense stiff BEM operator. Recent fractional and transient BEM developments [9,10,11,26,27,28,29] reinforce the relevance of boundary formulations, while the stiff predictor-corrector limitation identified in [6] explains why the implicit correction is timely.

Algorithm 2: Post-assembly operator stability certificate

Assemble the boundary-condition-eliminated DRBEM evolution matrix and compute its rightmost spectral abscissa. Estimate the departure from normality using the numerical range or selected resolvent norms when the operator is nonsymmetric or strongly nonnormal. Confirm that the dissipative spectrum lies inside the stability sector of the implicit fractional Adams–Moulton correction. Evaluate the nonlinear or explicitly treated Lipschitz scale and apply only the corresponding predictor restriction. If the resolvent margin deteriorates under refinement, refine the boundary interpolation, rebalance the dual-reciprocity point set, or reduce the explicit splitting before time integration. Record the spectral margin, conditioning indicator, and selected time-step policy as part of the reproducibility metadata.

Table 7 consolidates the principal verification outcomes into a compact publication-level evidence matrix. The purpose is to connect each numerical experiment to the analytical property it tests and to make the validation chain auditable without duplicating the detailed convergence tables.

9. Computational Complexity and Implementation Considerations

For $N_b$ boundary elements and $N_i$ DR interior points, the classical dense DRBEM assembly has the familiar $O\left(N_b^2\right)$ boundary integration cost plus the cost of forming and factorizing the RBF interpolation matrix. The Schur-complement operator $A_h$ is dense. For a fixed time step, the matrix $I-\omega A_h$ is factorized once; subsequent steps require $O\left(N_i^2\right)$ dense triangular solves.

Without history acceleration, fractional Adams accumulation costs vector operations and stores history data. SOE compression reduces the history dimension to recurrence states, leading to approximately history work and history storage [12]. For large-scale three-dimensional BEM, hierarchical matrices, adaptive cross approximation, fast multipole BEM, or iterative solvers can be used for the dense spatial operators; these accelerations are independent of the time-integration modification proposed here.

The method is especially attractive when: (i) the boundary is geometrically complicated relative to the volume, (ii) accurate boundary fluxes are important, (iii) remeshing a volume would be expensive, or (iv) the fractional model is coupled to boundary-dominated physics. It is less attractive when the governing operator has no convenient fundamental solution or when strongly heterogeneous coefficients make the dual-reciprocity interpolation excessively large.

10. Limitations and Extensions

Several limitations should be stated explicitly. First, the numerical implementation uses constant boundary elements and a simple global radial basis. Higher-order boundary interpolation and localized/compactly supported RBFs should improve conditioning and reduce the number of DR points. Second, the rigorous stability result is strongest for a time-independent linear stiff operator plus a Lipschitz explicit remainder. Fully nonlinear diffusion requires a nonlinear implicit treatment or a carefully designed frozen-Jacobian strategy. Third, direct history accumulation was used in the validation; SOE acceleration is formulated as an implementation extension and should be benchmarked separately in a long-time large-scale study. Fourth, corner singularities on nonconvex domains reduce spatial regularity, as visible in the pre-asymptotic behavior of Table 5. Singular elements, local refinement, or enriched boundary bases are natural improvements.

Three extensions are particularly promising. A multi-term or distributed-order Caputo model can be handled by summing history operators while keeping the same BEM spatial matrix. A coupled thermoelastic or bioheat system can be reduced blockwise by DRBEM and integrated with block-implicit Adams corrections. Finally, an adaptive stability controller can monitor the numerical range of the current linearized BEM operator and choose the explicit/implicit partition dynamically.

11. Conclusions

A new stability-centered boundary element development of the fractional Adams-Bashforth-Moulton methodology has been presented. DRBEM converts the Caputo FPDE into a dense matrix fractional ODE system in which the assembled spatial operator contains the geometry and boundary-condition information. Treating this operator implicitly in the fractional Adams-Moulton corrector removes the bounded negative-real-axis stability interval that restricts a fully explicit PECE approach, while retaining an Adams-Bashforth predictor for nonlinear and nonstiff terms.

The analysis provides a direct spectral stability certificate, a resolvent criterion for nonnormal/nonlinear cases, and a combined space-time error description. The numerical experiments are consistent with the theory: the disk benchmark gives near second-order DRBEM spatial convergence and second-order temporal convergence for smooth compatible data; the explicit PECE stability estimate becomes prohibitively small under BEM refinement; and the implicit formulation remains stable on a nonconvex L-shaped domain with a stiff left-half-plane spectrum.

The most important practical conclusion is that BEM and fractional Adams integration should not be coupled by simply applying an explicit predictor-corrector to the semidiscrete BEM equations. The stiffness of the boundary-element diffusion operator must be part of the design. The DRBEM-IABM formulation achieves this by making the spatial operator an explicit component of both the stability analysis and the implicit correction, providing a geometry-flexible and stiffness-robust platform for more complex fractional multiphysics models.