An error-certified matrix Mittag–Leffler perturbation method for weakly nonlinear fractional systems

1 Introduction

Fractional differential equations provide a natural description of systems whose current state depends on a continuum of past states. Such memory effects arise in viscoelasticity, anomalous transport, thermoelasticity, damping, electrical networks, biological dynamics, and many other areas [1–6]. For a constant matrix A, the matrix exponential e^*At* of classical dynamics is replaced by the matrix Mittag–Leffler propagator Eα,1(Atα). Stable scalar and matrix evaluations of Mittag–Leffler functions are therefore central to both analysis and computation [7–10].

We study weakly nonlinear Caputo systems

DCtαX(t)=AX(t)+G(t,X(t)), X(0)=X0, 0<α<1,(1)

where X(t) ∈ ℝ^*n*, A ∈ ℝ^*n*×*n*, and ε controls the nonlinear contribution. The exact variation-of-constants equation is

X(t)=Eα,1(Atα)X0+ε∫0t(t-s)α-1Eα,α(A(t-s)α)G(s,X(s))ds.(2)

Equation 2 suggests a perturbation expansion built around the exact fractional linear propagator rather than around an integer-order surrogate.

Perturbative and iterative approaches for fractional equations include the Adomian decomposition method (ADM), homotopy perturbation method (HPM), homotopy analysis method (HAM), variational iteration method (VIM), and fractional Picard iteration [11–16]. These techniques are useful, but they have different objectives and validation mechanisms. The contribution here is not the generic use of a power series in a small parameter. It is the combination of the exact matrix Mittag–Leffler propagator, a rigorously identified perturbation hierarchy, explicit truncation and residual-to-error estimates, singularity-resolved Volterra quadrature, and independent benchmarking against established fractional solvers.

Table 1 summarizes this positioning. The entries describe the typical implementation of each method; specialized variants may possess additional features.

Table 1. Positioning of the proposed method relative to common fractional perturbation and iteration frameworks.

The main advances are:

• (i) A general Fréchet-derivative formula for the perturbation coefficients;
• (ii) An explicit first-order error bound and a general analytic remainder estimate;
• (iii) An a posteriori estimate converting the Volterra residual into a solution-error bound;
• (iv) An exact transformation that removes the quadrature endpoint singularity;
• (v) Validation against refined Adams–Bashforth–Moulton predictor–corrector (ABM–PECE) solutions and an implicit Grünwald–Letnikov (GL) discretization;
• (vi) Three numerical examples, including perturbation-order, fractional-order, long-time, refinement, and computational-cost studies.

2 Preliminaries and Volterra formulation

2.1 Caputo derivative and matrix Mittag–Leffler functions

For 0 < α < 1, the Caputo derivative of a sufficiently regular scalar function is

DCtαx(t)=1Γ(1−α)∫0tx′(s)(t−s)αds.(3)

For vector-valued functions it is applied componentwise. The two-parameter Mittag–Leffler function is

Eα,β(z)=∑k=0∞zkΓ(αk+β), α>0,(4)

while, for a square matrix A,

Eα,β(Atα)=∑k=0∞AktαkΓ(αk+β).(5)

The linear nonhomogeneous Caputo system

DCtαX(t)=AX(t)+F(t), X(0)=X0,

has the variation-of-constants representation

X(t)=Eα,1(Atα)X0+∫0t(t-s)α-1Eα,α(A(t-s)α)F(s)ds.(6)

Define

Φα(t)=Eα,1(Atα), Kα(t)=tα-1Eα,α(Atα).(7)

Then Equation 1 is equivalent to the nonlinear Volterra equation

X(t)=Φα(t)X0+ε∫0tKα(t-s)G(s,X(s))ds.(8)

3 Matrix Mittag–Leffler perturbation hierarchy

3.1 Analytic expansion of the nonlinear term

We seek

X(t;ε)=∑m=0∞εmX[m](t).(9)

The following regularity assumption makes the expansion precise.

Assumption 3.1 (Analytic nonlinearity). For every t ∈ [0, T], the mapping X↦G(t, X) is analytic in a neighborhood of the unperturbed trajectory X^[0](t), and its Fréchet derivatives are continuous in t on that neighborhood.

Let

G(t,∑m=0∞εmX[m](t))=∑m=0∞εmF[m](t).(10)

The zeroth coefficient is

F[0](t)=G(t,X[0](t)),(11)

and, for m≥1,

F[m](t)=∑k=1m1k!∑j1+⋯+jk=mjℓ≥1DXkG(t,X[0](t))[X[j1](t),…,X[jk](t)].(12)

Thus, coefficient matching in Equation 8 yields

X[0](t)=Φα(t)X0,(13)

and

X[m+1](t)=∫0tKα(t-s)F[m](s)ds, m≥0.(14)

In particular,

X[1](t)=∫0tKα(t-s)G(s,X[0](s))ds,(15)

X[2](t)=∫0tKα(t-s)DXG(s,X[0](s))X[1](s)ds.(16)

The first-order approximation is

Xapp(t)=X[0](t)+εX[1](t).(17)

Remark 3.2. For α = 1, Φα(t)=eAt and Kα(t)=eAt, so Equation 14 reduces to the classical exponential-propagator perturbation hierarchy.

4 Existence, perturbation error, and residual certification

Let CT=C([0,T];ℝn) with norm

Define

KT=sup0≤t≤T∫0t||Kα(t-s)||ds.(18)

The weak singularity is integrable and KT<∞ on finite intervals.

Let X[0]=ΦαX0 and define the closed tube

BR={X∈CT:||X-X[0]||∞≤R}.(19)

Assume that on this tube

Theorem 4.1 (Contraction and unique solvability). If

then the Volterra operator

(TX)(t)=X[0](t)+ε∫0tKα(t-s)G(s,X(s))ds(23)

maps BR into itself and is a contraction. Consequently, Equation 8 has a unique solution in BR.

Proof. For X∈BR,

so T is a self-map. For X,Y∈BR,

Banach's fixed-point theorem completes the proof.

4.1 Explicit first-order truncation estimate

Theorem 4.2 (First-order perturbation error). *Under the assumptions of Theorem 4.2, the exact solution X and the first-order approximation (Equation 17) satisfy*

*Hence the first-order truncation error is rigorously O(ε²) for fixed* T.

Proof. From Equations 8, 15,

X-Xapp=ε∫0tKα(t-s)[G(s,X(s))-G(s,X[0](s))]ds.

Moreover,

Applying the Lipschitz estimate gives

The bound (Equation 24) is deliberately explicit and may be conservative because it uses global suprema over the full tube.

4.2 General perturbation-series remainder

**Theorem 4.3 (Analytic mth-order remainder)**. *Suppose the hypotheses of Theorem 4.1 hold uniformly for complex ε = ζ with |ζ| ≤ ρ, where ρLRKT<1, and suppose G(t, ·) is analytic. Then the solution X(·;ζ) is analytic in |ζ| < ρ, its Taylor coefficients are precisely those generated by Equation 14, and*

where

Mρ=sup|ζ|=ρ||X(·;ζ)||∞.

Proof. Uniform contractivity implies analytic dependence of the fixed point on the complex parameter. Cauchy's coefficient estimate gives ||X[k]||∞≤Mρρ-k. Summing the geometric tail from k = m+1 gives Equation 25.

This theorem distinguishes convergence of the perturbation series from convergence of a Picard sequence: the two expansions solve the same integral equation but are indexed by different quantities.

4.3 Residual-to-solution error estimate

For any continuous approximation Y, define its Volterra residual by

RV[Y](t)=Y(t)-X[0](t)-ε∫0tKα(t-s)G(s,Y(s))ds.(26)

Theorem 4.4 (A posteriori residual certification). Assume the exact solution X and an approximation Y belong to the Lipschitz tube and q < 1. Then

Proof. Subtracting Equation 26 from the exact equation gives

X-Y=ε∫0tKα(t-s)[G(s,X(s))-G(s,Y(s))]ds-RV[Y].

Therefore

and rearrangement yields Equation 27.

The residual is therefore not used as an unsupported proxy for accuracy. Under the contraction condition, it yields a rigorous error certificate; outside that regime it remains a consistency diagnostic and must be supplemented by an independent benchmark.

5 Singularity-resolved numerical implementation

5.1 Removal of the weak endpoint singularity

The original first correction contains

X[1](t)=∫0t(t-s)α-1Eα,α(A(t-s)α)G(s,X[0](s))ds.(28)

The substitution s = tr gives a factor (1−r)^α−1, which is integrable but singular at r = 1. Applying ordinary Gauss–Legendre quadrature directly to that integrand is not optimal. We instead define

y=(1-r)α, r=1-y1/α.(29)

Since

(1-r)α-1dr=-1αdy,

the singularity cancels exactly and

X[1](t)=tαα∫01Eα,α(Atαy)G(t(1-y1/α),X[0](t(1-y1/α)))dy.(30)

The integrand in Equation 30 is nonsingular on [0, 1], so a standard fixed Gauss–Legendre rule is justified. The same transformation is used in the residual integral.

5.2 Matrix-function evaluation and scaling

For the 2 × 2 oscillator studied below, distinct eigenvalues permit

Eα,β(Atα)=aβ(t)I+bβ(t)A.(31)

This identity is a convenience specific to the minimal polynomial of the example; it is not a restriction of the general method. For larger matrices one may use Schur–Parlett algorithms, contour methods, or Krylov approximations of matrix-function actions [7, 9, 17].

The scalar coefficients are precomputed on a fine grid and represented by cubic splines. All numerical experiments use a 64-point Gauss–Legendre rule unless otherwise stated. The source code evaluates the scalar Mittag–Leffler series until the newest term is below a relative tolerance and raises an exception if convergence is not attained.

5.3 Independent benchmark solvers

The ABM–PECE benchmark follows the fractional Adams method of Diethelm, Ford, and Freed [18, 19]. For

DCtαY=f(t,Y), Y(0)=Y0,

the predictor is

Yn+1P=Y0+hαΓ(α+1)∑j=0nbj,n+1f(tj,Yj),(32)

where

bj,n+1=(n+1-j)α-(n-j)α.

The corrector is

Yn+1=Y0+hαΓ(α+2)[f(tn+1,Yn+1P)+∑j=0naj,n+1f(tj,Yj)],(33)

with the standard fractional Adams weights. A refined solution with h = 0.00125 or 0.0025 is used as ground truth, depending on the experiment.

For a second independent comparison, an implicit GL discretization is used:

h-α∑k=0nωk(Yn-k-Y0)=f(tn,Yn), ω0=1, ωk=(1-α+1k)ωk-1,(34)

with the nonlinear algebraic system solved at every step [20, 21].

The reported error indicators are

E∞=maxj||Xapp(tj)-Xref(tj)||,(35)

ERMS=(1N+1∑j=0N||Xapp(tj)-Xref(tj)||2)1/2,(36)

Erel=E∞1+maxj||Xref(tj)||.(37)

The complete Python reproduction script, generated numerical tables, machine-readable data, and publication-quality figures are provided in the Supplementary Material.

6 Example 1: exact manufactured verification

Consider

xDCtα(t)=−x(t)+εgα(t), x(0)=1,(38)

where

gα(t)=t2+2t2-αΓ(3-α).(39)

The exact solution is

x(t)=Eα,1(-tα)+εt2.(40)

Because the forcing is independent of x, the perturbation expansion terminates after the first correction. Thus this problem isolates errors caused by Mittag–Leffler interpolation and Volterra quadrature.

For α = 0.80, ε = 0.20, and T = 3, Table 2 gives the errors. The maximum actual error is 6.025 × 10⁻⁷, while the residual is at roundoff level because the same accurately transformed integral operator is evaluated on the numerical approximation. Figure 1 compares the exact solution with the first-order matrix Mittag-Leffler approximation; the two curves are visually indistinguishable at the displayed scale.

Table 2. Exact manufactured verification.

7 Example 2: fractional Duffing-type state-space system

7.1 Problem and explicit propagator

We consider

DCtαX(t)=AX(t)+εG(X(t)), X(0)=X0,(41)

with

A=(01-1-0.30), X0=(10), G(X)=(0-x13).(42)

This is a fractional Duffing-type state-space model. For 0 < α < 1, it should not automatically be rewritten as a scalar equation involving a sequential derivative of order 2α, because Caputo derivatives do not generally satisfy a semigroup property.

The eigenvalues of A are

λ±=-0.15±0.988686i.

For distinct eigenvalues,

bβ(t)=Eα,β(λ+tα)-Eα,β(λ-tα)λ+-λ-,(43)

aβ(t)=λ+Eα,β(λ-tα)-λ-Eα,β(λ+tα)λ+-λ-,(44)

and

Eα,β(Atα)=(aβ(t)bβ(t)-bβ(t)aβ(t)-0.30bβ(t)).(45)

The unperturbed solution is

X[0](t)=(a1(t)-b1(t)),(46)

and the first correction is

X[1](t)=-∫0t(t-s)α-1(bα(t-s)aα(t-s)-0.30bα(t-s))a13(s)ds.(47)

7.2 Base-case validation against PECE

The base parameters are

α=0.85, ε=0.10, T=3.(48)

The approximation is sampled at h = 0.01, while the reference solution uses ABM–PECE with h = 0.0025. Table 3 reports the actual error and the independently evaluated Volterra residual. Figures 2–4 compare the state components, the phase curves, and the pointwise actual error with the Volterra residual, respectively.

Table 3. Duffing base-case error and residual validation.

The actual maximum error is 2.131 × 10⁻³ and the residual is 2.176 × 10⁻³. On the tube |x₁| ≤ 1, one may take L_*R* = 3 and M_*R* = 1. Numerical evaluation gives

KT=2.13653, q=εLRKT=0.64096<1.

The a posteriori theorem therefore yields

which safely contains the measured error. The explicit a priori bound (Equation 24) is 1.369 × 10⁻¹ and is much more conservative, as expected from global tube constants.

7.3 Perturbation-order study

To test the claimed order directly, the same X^[0] and X^[1] are used for

ε∈{0.025,0.05,0.10,0.20}.

For consecutive perturbation values, the observed order is

pi=log(E∞,i/E∞,i-1)log(εi/εi-1).(50)

Table 4 and Figure 5 show orders 1.991, 1.960, and 1.918. The gradual reduction at ε = 0.20 is consistent with the increased importance of higher powers.

Table 4. Perturbation-parameter sensitivity and observed first-order truncation rate.

7.4 Fractional-order sensitivity

Table 5 reports results for

α∈{0.65,0.75,0.85,0.95}.

remain of order 10⁻³. The error increases moderately as α approaches one for this parameter set, but no loss of stability is observed. Figure 6 displays the corresponding dependence of the actual error and Volterra residual on the fractional order.

Table 5. Sensitivity to the fractional order for ε = 0.10 and T = 3.

7.5 Long-time behavior and contraction margin

The interval is extended to T = 6 and T = 10. Table 6 shows that the maximum error remains controlled, while KT and q increase. At T = 10, q = 0.869, so the contraction-based certificate remains valid but has a narrower margin. This behavior illustrates an important limitation: even when the observed trajectory remains accurate, a global contraction estimate can become ineffective as the memory horizon grows.

Table 6. Long-time Duffing study.

CPU is the single-threaded perturbation cost, including propagator precomputation.

7.6 Reference refinement and method comparison

A PECE refinement study is included to verify that the selected ground truth is numerically reliable. The observed orders approach 1+α = 1.85, in agreement with fractional Adams theory [19].

Table 7 reports the ABM-PECE refinement errors, observed orders, and CPU times, while Table 8 compares the first-order matrix Mittag-Leffler method, ABM-PECE, and the implicit Grünwald-Letnikov method at the nominal output step h = 0.01. The perturbation approximation is not intended to replace a high-accuracy time-stepping solver when a single trajectory is required: PECE is more accurate and slightly faster in this small test. Its advantage is structural and parametric: once X^[0] and X^[1] are computed, all first-order predictions for new ε values require only the affine combination X^[0]+εX^[1]. By contrast, PECE and GL must repeat the complete memory calculation.

Table 7. ABM–PECE refinement against a solution with h = 0.00125.

Table 8. Method comparison for the Duffing base case.

Timings are single-threaded Python 3.13.5 measurements in the same execution environment.

8 Example 3: nonlinear fractional logistic equation

To demonstrate that the procedure is not restricted to the oscillator, consider

DCtαx(t)=−0.40x(t)+εx(t)(1−x(t)), x(0)=0.70,(51)

with

α=0.80, ε=0.25, T=5.

The nonlinear term is analytic and the first correction is obtained from Equation 30 with scalar propagator Eα,β(-0.40tα). A PECE solution with h = 0.0025 is used as reference. Table 9 reports the actual error and Volterra residual obtained for the nonlinear fractional logistic example.

Table 9. Error and residual for the nonlinear fractional logistic example.

The maximum relative error is 9.51 × 10⁻³, below one percent under the adopted normalization. The larger absolute error compared with the Duffing example is expected because ε = 0.25 and the nonlinear term acts throughout the five-unit interval. Figure 7 compares the first-order matrix Mittag-Leffler approximation with the refined PECE reference solution.

9 Computational complexity, applicability, and limitations

For N output times and a fixed q_*G*-point quadrature, a direct first-order implementation costs approximately O(Nq_*G*) matrix-function actions after propagator precomputation. For a dense n×n matrix this action can cost O(n³) if a full matrix function is formed, while Krylov methods reduce the cost when only its action on a vector is required. Classical PECE and GL history methods have a direct O(N²) memory-convolution cost unless fast convolution or short-memory acceleration is introduced [22–24].

The matrix perturbation approach is most useful when:

• (i) A captures the dominant dynamics;
• (ii) ε is sufficiently small for a low-order truncation to be meaningful;
• (iii) The propagator can be precomputed or efficiently applied;
• (iv) Many values of ε, initial data, or weak nonlinear coefficients are explored with the same linear operator.

Its limitations are equally important. Accuracy may deteriorate for strong nonlinearities, large Lipschitz constants, unstable linear matrices, long time intervals, or trajectories leaving the tube used in the estimates. The a priori bound can be conservative, and the residual certificate requires q < 1. For large matrices, the 2 × 2 interpolation shortcut used in the Duffing example must be replaced by robust matrix-function algorithms. Higher-order terms also become algebraically complex because of the multilinear derivatives in Equation 12. These points delimit the method's applicability and prevent the residual from being interpreted as a universal substitute for a reference solution.

The numerical results also show that residual and actual error need not coincide exactly. In the Duffing case they are close, while in the logistic example the actual error is somewhat larger. The residual theorem explains their relationship only after multiplication by the stability factor (1−q)⁻¹. Consequently, the recommended validation protocol is: residual evaluation, perturbation-order verification, and at least one refined independent numerical comparison.

10 Conclusion

A matrix Mittag–Leffler perturbation framework has been developed for weakly nonlinear Caputo systems. The proposed framework addresses the principal theoretical and numerical gaps of a residual-only approach. The perturbation coefficients are defined by an analytic Fréchet hierarchy; the first-order approximation has an explicit O(ε²) error estimate; a general analytic-series remainder is provided; and the Volterra residual yields an a posteriori solution-error bound when the contraction factor is below one.

The weak singularity of the memory integral is removed exactly by a power substitution, allowing standard Gauss–Legendre quadrature on a nonsingular integrand. An exact manufactured problem verifies the implementation to approximately 6 × 10⁻⁷. For the Duffing state-space system, the actual error and residual are both approximately 2.1 × 10⁻³, and the observed perturbation orders are close to two. PECE refinement, GL comparison, fractional-order sensitivity, long-time tests, and CPU data provide independent evidence of accuracy and clarify that the principal advantage is reuse of the exact linear fractional propagator in weakly nonlinear parameter studies. A nonlinear logistic example confirms applicability beyond oscillators.

Future work will address second- and higher-order automatic differentiation, fast history compression, adaptive interpolation of matrix-function actions, and Krylov implementations for large sparse systems.

Generative AI statement

The author(s) declared that Generative AI was used in the creation of this manuscript. During manuscript preparation and the production-correction stage, the authors used ChatGPT (OpenAI) for language refinement, organization of the structured abstract, and preparation of responses to the production queries. The authors critically reviewed, edited, and approved all AI-assisted text and take full responsibility for the final article. The mathematical derivations, numerical computations, data, figures, and scientific conclusions were generated and verified by the authors and the supplied reproducibility code.

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

Publisher's note

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

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fams.2026.1899674/full#supplementary-material

Python source code, generated numerical data, CSV and LaTeX tables, figures in PDF and PNG formats, and a machine-readable JSON summary of the reported metrics. Supplementary material at the source

Footnotes

Edited by: Khalid Hattaf, Centre Régional des Métiers de l'Education et de la Formation (CRMEF), Morocco

Reviewed by: Emad Awad, Alexandria University, Egypt

Livija Cveticanin, University of Novi Sad, Serbia

Muhammad Umer, University of Sargodha, Pakistan