Computer Science and Mathematics

Abstract: Exposing foams stabilized by photoswitchable surfactants to UV light induces changes in surface surfactant concentration, leading to significant alterations in foam behaviour such as the generation of Marangoni flow and change in foam drainage patterns. The occurrence of Marangoni flow can be observed when either all elements of the foam or only their films are exposed to UV light. Conversely, changes in foam drainage occur when a macroscale portion of a foam column is exposed to UV light. To explore these phenomena, numerical models are developed and validated using experimental data. These models simulate the scale and profile of Marangoni flow from foam networks to films as well as the drainage flow within the foam network. Microscale findings demonstrate that Marangoni flow can be controlled by adjusting the intensity and duration of UV light exposure. Macroscopically, the drainage profile in exposed foam regions undergoes significant changes with varying UV intensity. Furthermore, beyond a certain threshold, the foam drainage reverses direction, contrary to gravity. The effect of foam interfacial mobility on the reversed drainage of both interior and exterior foams is analyzed. The findings provide a potential tool to control foam drainage behaviour without the need to modify other variables.

Abstract: Accurate and computationally efficient prediction of radioactive plume dispersion is essential for rapid risk assessment and emergency response in the event of accidental atmospheric releases. This study presents a three-dimensional, mass-conservative numerical framework for simulating mesoscale transport and transformation of radionuclides under varying meteorological conditions. The proposed model is based on the advection–diffusion equation and incorporates key physical processes, including turbulent diffusion, gravitational settling, precipitation washout, radioactive decay, surface absorption, and re-emission. A hybrid numerical scheme is developed that combines a semi-implicit Crank–Nicolson discretization for diffusion terms with a second-order Total Variation Diminishing (TVD) scheme employing a Van Leer limiter for advection. To ensure computational efficiency and stability, the governing equations are solved using an alternating direction implicit (ADI) strategy, leading to tridiagonal systems efficiently handled via the Thomas algorithm. In addition, a modified parameterization of turbulent diffusion coefficients is introduced, extending classical Pasquill–Gifford formulations to better represent mesoscale atmospheric behavior with bounded dispersion characteristics. Computational experiments simulating iodine-131 release under different atmospheric stability classes and meteorological conditions demonstrate that the proposed framework reproduces physically consistent plume structures and maintains strict mass conservation with negligible numerical errors. The results highlight the dominant role of advective transport and precipitation scavenging under moderate wind conditions, as well as strong localization effects under stable atmospheric stratification. The proposed model provides a practical compromise between the simplicity of Gaussian plume approaches and the high computational cost of full-physics atmospheric models, making it suitable for rapid engineering calculations and integration into decision-support systems for radiological emergency management.

Abstract: We develop a unified and fully audited analytic hierarchy for the restricted weighted Goldbach sum$$R_{a,q}(N) := \sum_{\substack{p_1+p_2=N \ p_1 \equiv a \pmod{q}}} (\log p_1)(\log p_2), \qquad q \geq 1,\ \gcd(a,q)=1,$$with expected main term $M_{a,q}(N) := C_2, S(N), N/\varphi(q)$, and exceptional set $E_{a,q}(X) := {N \leq X,\ N\ \text{even} : R_{a,q}(N) = 0}$. The paper consolidates and supersedes preprint version 3, integrating results from Papers 1, 9 and 14 of the Anderson Series, with all documented corrections applied.The unconditional core establishes three nested levels. Level 1 is an effective almost-all theorem via the standard $L^4$ minor-arc route, with explicit constant $K = 2C(1,4) \leq 38.82$. Level 1.5 is a sub-exponential exceptional-set bound $\#\mathcal{E}_{a,q}(X) \ll_{q} X \exp\left(-\frac{\sqrt{\log X}}{R}\right)$ with Stechkin's constant $R = 9.6459$, proved unconditionally by absorbing any potential Siegel zero into a modified main term. Level 1.5+ is a Hölder minor-arc refinement giving the improved constant $K_{\text{new}} \leq 9.80$ and, for moduli $q \leq 200$ certified free of Siegel zeros, an unconditional pointwise sub-exponential bound with $C(4) \leq 120$ and $\log N_0(4) \leq 42$.Three structural obstructions (Double-Pole, Borel–Cantelli, ETK Dimensional Explosion) formally retract three classical routes to unconditional finiteness. Under DH and GRH, conditional hierarchies (with $\theta(A) = 1 - 2/(A+2)$ and $\log N_0(4) = 45.93$) are recorded. The Gowers–Spectral Bridge gives conditional finiteness under the Uniform Spectral Gap (USG) hypothesis with effective threshold $N_0(4) \leq 10^{16}$.The open sub-lemma (Proposition 13.5) connecting USG to the Montgomery Pair Correlation Conjecture has been resolved conditionally. The complete logical chain $\text{Strong Montgomery-GUE} \Rightarrow \text{USG} \Rightarrow R_{a,q}(N) > 0$ for all $N \geq N_0$ is now fully established (conditional on Strong Montgomery-GUE, which is strictly stronger than the standard weak form of Conjecture 13.10). See Sections 13.5–13.6.[HONEST CAVEAT] The proof of Theorem 13.15 invokes EAC (Controlled Additive Energy) as an intermediate step. By Theorem 13.16, EAC is equivalent to Strong Montgomery-GUE, which is strictly stronger than the standard Montgomery Pair Correlation Conjecture (Conjecture 13.10). The logical chain should therefore be read as: Strong Montgomery-GUE $\Rightarrow$ USG $(c=1)$ $\Rightarrow$ finiteness of $E_{a,q}$. The standard weak Montgomery conjecture alone does not suffice; see Section 13.6, Obstacle 1.

Abstract: Scientific inquiry relies on causal explanation, whereby phenomena are understood through the antecedent conditions, interactions and laws producing the current configuration of a system. We introduce a complementary approach termed Future Compatibility (FC), which evaluates a present state not only through its causal history, but also through the set of future states structurally accessible from it. The central premise is that scientific understanding may be enriched by characterizing the opportunities and constraints associated with a given state. Rather than focusing exclusively on origins or on the prediction of a specific outcome, FC quantifies the repertoire of coherent futures available from a present configuration and their stability under perturbation. To explore this approach, we implemented simulations of evolving networked systems undergoing branching developmental pathways and external perturbations. For each system state, a set of novel observables describing future accessibility, entropy, robustness and persistence was computed from admissible trajectories satisfying predefined coherence constraints. These compatibility landscapes exposed organizational properties that were not captured by causal reconstruction alone, providing additional information regarding adaptability, resilience and developmental potential. Potential applications include developmental biology, evolutionary theory, ecology, neuroscience, complex systems research and the assessment of resilience and adaptability in technological, organizational and socioeconomic networks.

Abstract: We develop P versus NP as the feasibility instance of the horizon principle over a finite relational substrate. The algebraic carrier is complete as a totality, yet a bounded internal observer has a finite comprehension horizon. The incompleteness migrates from the undecidability of truth to that horizon inaccessibility. Computation is counting along a representation: the hyperoperation ladder (succession, addition, multiplication, exponentiation) carried by the substrate's four cardinal representation charts. A certificate is a near-representation placing a remote target within the horizon: NP asks that one exists, P that it be forward-found, and P=NP iff finding is as feasible as checking. Geometrically the carrier is a shell with the observer at its pole and a computation a geodesic, so P=NP asks whether a geodesic that provably exists is forward-findable. The reading reproduces, under exact substitutions, proof complexity, automatizability, and one-way functions, and grounds the find-check asymmetry in the substrate's scale/discrete-logarithm map and single time arrow. P versus NP is the computational face of the horizon clause: its uniform certificate is an Ω-hard residue, decided by the complete totality yet below the bounded observer's horizon. One-way functions yield the separation; the descent's one-wayness is that residue. The reading is consistent with the relativization, natural-proofs, and algebrization barriers. Every exact claim is verified in finite-field or cyclotomic arithmetic.

Abstract: We introduce a generalized fixed-point framework based on $d$-$\psi$-$\varepsilon$-contractions and weak $d$-$\psi$-$\varepsilon$-contractions, where the contraction metric and the metric in which the space is complete may be different and the comparison function belongs to the class $\Psi_0(\varepsilon)$. This setting extends several classical fixed-point principles and yields existence, uniqueness, localization, and convergence results for Picard iterations. The theory is applied to nonlinear integral equations and to their quadrature approximations. By introducing suitable invariant sets and a \(\psi\)-\(\varepsilon\)-max inequality, existence and uniqueness results are obtained under assumptions that are substantially different from classical Lipschitz-type conditions. The developed framework is further extended to quadrature integral equations generated by numerical integration formulas. Sufficient conditions are established for the existence and uniqueness of solutions as well as for the convergence of the associated Picard sequences. The theoretical results guarantee convergence of the discrete iterations under appropriate assumptions, thereby providing a direct connection between fixed-point theory and numerical computation. Several examples, including Chandrasekhar-type integral equations and nonlinear weighted integral equations with singular kernels, are presented to illustrate the applicability and effectiveness of the proposed approach. Numerical experiments confirm the theoretical findings and demonstrate the accuracy of the resulting quadrature-Picard schemes.

Abstract: Fracture of ceramic femoral heads in total hip arthroplasty is a rare but catastrophic complication requiring urgent revision surgery. Most finite element studies are limited to static loading and do not capture the dynamic behavior of ceramic components under impact conditions generated during stumbling or falling. In the present study, a parametric explicit dynamic analysis was performed using LS-DYNA (version 960) to determine the impact fracture thresholds of alumina (Al₂O₃) and yttria-stabilized zirconia (ZrO₂, Y-TZP) femoral heads. An axisymmetric finite element model of a 32 mm ceramic femoral head articulating with a ceramic liner within a Ti-6Al-4V acetabular shell was developed. Ceramic behavior was described using the Johnson–Holmquist JH-2 damage constitutive model. The viscoelastic bone stock response was represented by a Winkler foundation of discrete spring-dashpot elements (stiffness 50–500 N/mm, damping 0–1.0 N·ms/mm). Impact velocity was varied from 0.01 to 0.45 mm/ms, consistent with velocities recorded by instrumented implant telemetry during stumbling. Fracture was identified by three concurrent criteria: effective plastic strain, maximum principal stress, and inflection of the internal energy–time curve. For Al₂O₃, the critical fracture velocity was 0.08 mm/ms under rigid fixation and 0.05 mm/ms with a viscoelastic foundation. The ZrO₂ femoral head did not fracture at any velocity tested; at V ≥ 0.20 mm/ms, the Ti-6Al-4V neck underwent plastic deformation as a competing failure mode while the ceramic head remained intact. Foundation stiffness and damping had no influence on fracture outcome across the clinically relevant range, indicating inertia-dominated fracture mechanics. These results provide quantitative fracture thresholds to support comparative evaluation of ceramic hip implant components.

Abstract: We present a self-contained, fully audited analytic hierarchy for the restricted weighted Goldbach sum$$R_{a,q}(N) ≔ \sum_{\substack{p_{1} + p_{2} = N \\ p_{1} \equiv a\ (mod\ q)}}^{}(\log p_{1})(\log p_{2}),\quad\quad q \geq 1,\mspace{6mu} \gcd(a,q) = 1,$$with expected main term $M_{a,q}(N): = C_{2}S(N)N/\varphi(q)$, and exceptional set $\mathcal{E}_{a,q}(X): = \{ N \leq X,\, N\text{ even}:R_{a,q}(N) = 0\}$.Part I (Unconditional core). We establish, with certified numerical constants, an effective almost-all theorem with explicit threshold $K = 2C(1,4) \leq 38.82$; a uniform minor-arc $L^{4}$ bound with $\kappa_{safe} \leq 4.40$; the exact diagonal second-moment constant $G/(2\varphi(q))$; and the structural rigidity results (gap bound with exponent $0.525$, non-consecutiveness, three-term progression avoidance, and $\parallel 1_{\mathcal{E}_{a,q}} \parallel_{U^{2}\lbrack 1,X\rbrack} \rightarrow 0$).At the Level 1.5 milestone, we embed the full proof of the sub-exponential exceptional-set bound$$\#\mathcal{E}_{a,q}(X) \ll_{q}X\exp\left( - \sqrt{\frac{\log X}{R}} \right),\quad R = 9.6459\text{ (Stechkin)},$$proved unconditionally by absorbing any Siegel zero into the modified main term $M_{a,q}^{mod}(N)$.We additionally record the convergence and Cesàro identification of the amplification factor $$S_{\infty} ≔ \prod_{\mathcal{l} > 2,\mathcal{\, l \in P}}^{}\left( 1 + \frac{1}{\left( \mathcal{l} - 1 \right)\left( \mathcal{l} - 2 \right)} \right) = 1.74272535539183\ldots,$$which governs the shifted-prime subsequence $\{ p + 1:p\text{ prime}\}$ and worsens the explicit threshold constant to $K_{new} \approx 51.3 = \sqrt{S_{\infty}}\, K$, a phenomenon we term amplification penalty.Part II (Structural obstructions). Three classical routes to unconditional finiteness are formally refuted. The Double-Pole Convolution Obstruction shows that the binary problem is governed by $( - L'/L)^{2}$, not $- L'/L$, so a fixed Siegel zero $\beta_{1} < 1$ contributes only $O(N^{2\beta_{1} - 1}) = o(N)$. The Borel–Cantelli Divergence Barrier proves $\mu(A_{N}) \gg 1/N$ under the Linear Independence Conjecture, so $\sum_{N}^{}\mu(A_{N}) = \infty$. The ETK Dimensional Explosion shows the Erdős–Turán–Koksma error factor $3^{k(N)} \rightarrow \infty$ defeats Baker-type bounds. Part III (The Gowers–Spectral Bridge). Under the Uniform Spectral Gap hypothesis (USG), a statistical regularity condition on the zero-sum graph of Dirichlet $L$-functions strictly weaker than GRH, the exceptional set $\mathcal{E}_{a,q}$ is finite with an effectively computable threshold $N_{0}(q)$. For $q = 4$, one obtains $N_{0}(4) \leq 10^{16}$ under USG with effective phase dimension $d_{eff} = 4$.Every statement carries one of the epistemic labels \[PROVED\], \[COND. PROVED, $H$\], \[OPEN\], \[RETRACTED\], or \[HONEST CAVEAT\].

Abstract: In the single-depot multi-traveling salesman problem, traditional depot location methods often overlook task balance among traveling salesmen, leading to excessive load on certain units and compromising overall operational efficiency. To address this issue, this paper proposes an optimized depot location method based on clustering and multi-task balancing. The core contribution lies in the design of a multi-weight adaptive depot optimization method. This approach clusters city nodes into multiple groups through cluster analysis and dynamically synthesizes direction vectors using information such as the number of samples within each cluster and the convex perimeter. It iteratively optimizes depot locations, minimizing the total path length while enhancing load balance across all traveling salesman routes. Additionally, a “divide-and-conquer” strategy decomposes the complex MTSP into multiple parallel TSP subproblems, which are then efficiently solved using Or-Tools. A comprehensive evaluation framework is introduced, incorporating Total-Sum distance, Min-Max distance, Distance-balancedness, Cluster separability, Robustness, and Running time. Experimental results on the TSPLIB standard dataset demonstrate that the proposed method exhibits significant advantages over various traditional clustering algorithms in both route optimization and route balancing, validating its effectiveness and practicality. The method's robust performance provides a reliable solution for real-world applications such as logistics distribution, further highlighting its practical value.

Abstract:

Nonlinear matrix equation X^p - ^m∑_i=1A^T_iX^-1A_i=Q has extensive applications in control theory, ladder networks, dynamic programming and stochastic filtering. In this paper, we propose a weighted average iterative algorithm to solve this nonlinear matrix equation. Based on the basic characteristics of the Thompson metric space, the convergence and error estimation formula of the algorithm are proved. Numerical experiments to illustrate the feasibility and effectiveness of the proposed algorithm are given.

Abstract: This study presents a comparative numerical investigation of second-order and fourth-order Runge–Kutta methods for solving chaotic dynamical systems. The Lorenz, Genesio–Tesi, and Rossler systems are considered because of their nonlinear behavior and high sensitivity to initial conditions. The numerical schemes investigated include the Midpoint, Improved Euler, Ralston, and fourth-order Runge–Kutta (RK4) methods. The performance of the methods is evaluated in terms of convergence behavior, numerical accuracy, stability characteristics, and computational cost. Stability analysis of each chaotic system is carried out through equilibrium point determination and Jacobian eigenvalue analysis. Numerical simulations are implemented in MATLAB and comparisons are performed using different step sizes. The results indicate that all numerical methods converge as the step size decreases; however, the RK4 method consistently provides significantly smaller errors and improved stability properties compared with the second-order schemes. The findings further demonstrate that higher-order numerical integration methods provide superior performance for highly sensitive chaotic systems where accuracy and reliability are essential.

Abstract: Physical infrastructure failures expose a structural blind spot in conventional financial risk systems: when the electricity supply becomes binding, markets enter metabolic arrest, a state where information-processing capacity collapses, and standard correlation-based metrics misread rising cross-asset correlations as stability. This paper introduces TENSORnet (Temporal Entropy-aware Network for Systemic Onset Recognition), a physics-informed computational protocol that fuses a VIX-calibrated inhibitory gate γL(t) (NLS-fitted, R2 = 0.30) with cross-asset Shannon entropy to detect infrastructure-induced stress before it materialises in price-based indicators. Applied to a temporal cross-asset graph of 2,838 JSE trading days (05 January 2015 to 29 April 2026) across seven asset classes under South Africa’s load-shedding crisis, TENSORnet achieves Precision = 100.0%, Recall = 85.8%, F1 = 92.4%, AUC = 1.000, and zero false positive alarms (Stage 3+ definition) on 1,830 out- of-sample days, with a mean lead time of 17 calendar days (408 hours) before fragile regime onset, outperforming all benchmarks, including XGBoost (F1 = 71.3%). Ablation confirms the physics-informed gate as the dominant architectural component (ΔF1 = −92.4 pp on removal): statistical learning alone recovers nothing (AUC = 0.469). The Densification Paradox, rising cross-asset correlation with falling entropy under stress (r = −0.468, p < 0.001), is confirmed empirically for the first time in cross-asset data. A Thermodynamic Port-Hamiltonian Neural ODE (TPH-NODE) extension grounds metabolic arrest in the second law of thermodynamics.

Abstract: Differential Causal Networks (DCNs) were introduced to represent changes between two causal networks inferred under different conditions. In their original use, however, DCNs remain pairwise objects: each differential graph summarizes rewiring within a single system, while common differential structures shared across many systems remain implicit. We introduce a methodological framework for the local alignment of DCNs aimed at detecting recurrent rewiring motifs, that is, small directed differential subnetworks that reappear across multiple systems under the same contrast. The proposed framework transforms each system-specific comparison into a signed directed differential graph, preserves both edge direction and change type, and searches for approximate local correspondences rather than enforcing a full-network mapping. The method consists of four steps: construction of signed DCNs, extraction of differential seeds, pairwise local alignment by seed-and-extend, and progressive multiple alignment to build consensus motifs. We define a score that combines node compatibility, differential-edge conservation, directional consistency, and recurrence support, and we complement the alignment procedure with null-model testing and robustness analysis. The result is a collection of consensus local differential modules ranked by recurrence, confidence, and statistical significance. In this formulation, DCNs become comparable units in a higher-order analysis whose goal is not merely to describe pairwise causal change, but to identify the same local rewiring logic reused across multiple systems.

Abstract: Optimisation algorithms play an important role in the solution of nonlinear engineering design problems, particularly where objective functions exhibit complex, nonconvex, and potentially multimodal behaviour. Classical gradient-based methods, including the Method of Moving Asymptotes (MMA) and Sequential Quadratic Programming (SQP), are widely recognised for their computational efficiency and rapid local convergence; however, their performance may be sensitive to the presence of local extrema. In contrast, metaheuristic approaches such as Particle Swarm Optimisation (PSO) generally provide enhanced global exploration capabilities, albeit often at significantly greater computational expense. This study presents a preliminary investigation of a hybrid optimisation framework termed the Constrained Adaptive Model-based Exploration Optimiser (CAMEO). The proposed approach combines bounded stochastic exploration with constrained local refinement in an attempt to improve robustness within multimodal optimisation landscapes whilst retaining the efficiency associated with deterministic optimisation methods. The performance of the proposed framework was examined using a series of benchmark optimisation problems and compared against MMA, SQP, and PSO. The numerical results indicate that CAMEO is capable of attaining solutions closer to the global optimum in several test cases, whilst maintaining stable convergence characteristics.

Abstract: Modern machine learning systems can achieve remarkable predictive performance. Nevertheless, in several fields, this is not enough to produce acceptable solutions as we need formal guarantees of robustness, fairness, and interpretability. Most existing approaches treat these properties separately or introduce them through external constraints, which makes their interaction difficult to analyze. In this work, we develop a unified variational perspective that incorporates these requirements directly into the learning objective. Concretely, we model learning as the minimization of a composite functional that combines predictive risk, regularization, and additional terms that capture robustness, fairness, and interpretability. This viewpoint allows us to study these properties within a single mathematical framework. Under standard assumptions, we prove the existence of minimizers and show that the resulting solutions are Pareto-optimal for the associated multi-objective problem. We illustrate the framework using examples based on adversarial and distributional robustness, statistical fairness criteria, and a notion of interpretability. The analysis points out the trade-offs that inevitably arise. We also examine statistical aspects of the proposed objective and show that classical generalization guarantees can still be obtained under appropriate conditions. The resulting framework provides a flexible basis for designing reliable learning systems.

Abstract: We present and analyze a new weighted family of iterative methods for solving systems of nonlinear equations. The proposed schemes are constructed as a generalization of the fifth-order method of Singh et al. by incorporating appropriate weight functions into the correction step, thereby generating a flexible class of methods that includes the original scheme as a special case. Sufficient conditions on the weight functions are established to guarantee fifth-order local convergence. Several admissible choices are presented to illustrate the versatility of the family. The practical performance of the proposed variants is investigated on a collection of large-scale nonlinear systems. Furthermore, the family is applied to the nonlinear algebraic system obtained from the finite-difference discretization of a stationary one-dimensional viscous Burgers problem. Numerical experiments indicate that the proposed methods provide a competitive and accurate alternative for solving nonlinear systems of this type.

Abstract: This paper develops a robust numerical scheme based on a frame collocation method for solving multi-term fractional ordinary differential equations (FODEs) whose solutions exhibit multiple singularities at the origin. To adaptively capture the singular behavior, we construct a hybrid basis function frame by combining shifted fractional Legendre polynomials. An efficient computational formula for the Caputo fractional derivative is derived, which transforms the original problem into a nonlinear algebraic system at the collocation points. The resulting system matrix is severely ill-conditioned due to the redundancy of frame, to mitigate this issue, we employ truncated singular value decomposition (TSVD) regularization, thereby enabling stable and high-precision solutions. Extensive numerical experiments on several benchmark problems, including the fractional Bagley–Torvik equation, linear multi-term FODEs, and nonlinear cases, demonstrate that the proposed method achieves exponential convergence rates. Notably, when the singular exponent of the solution matches a tunable parameter $\delta$ in the basis functions, superconvergence is observed, significantly outperforming standard spectral methods. Compared with traditional spectral approaches, the proposed frame collocation framework retains spectral accuracy while exhibiting superior capability in handling complex singular structures, providing a powerful and reliable tool for high-precision simulations of multi-term fractional differential equations.

Abstract: APD is widely adopted in the management of end-stage renal disease (ESRD) and offers flexi-bility and improved quality of life, but bacterial infections, particularly peritonitis, are still a major constraint, which frequently results in hospitalization, catheter failure, and hemodialysis. Early diagnosis is important but difficult because of the non-specific clinical manifestations and delays related to the traditional diagnostic techniques like culture-based analysis. “To overcome these restrictions, this paper suggests a new explainable machine learning model to early identify bacterial infections in APD patients based on multimodal data streams, such as clinical, lab, and time-series dialysis data, to identify both fixed and dynamic infection onset patterns”. The framework uses a hybrid characteristic of feature engineering, which is a combination of statistical selection techniques and clinically relevant indicators to improve predictive performance, and Supervised learning models of high accuracy like the Random Forest, SVM, and Gradient Boosting are applied. One of the contributions of this work is the incorporation of explainable artificial intelligence through SHAP that leads to a clear interpretation of model predictions and the determination of key risk factors that will affect the development of the infection and thus enhance clinical trust and usability. The experimental findings indicate that the given approach greatly enhances the accuracy of early detection as compared to the conventional ones, allowing timely intervention, minimizing complications, and improving the overall outcomes of the treatment, which underscores its potential as a scalable and clinically applicable decision support system to manage APD.

Abstract: We develop a fast Chebyshev spectral collocation method for a coupled system of nonlinear Klein–Gordon equations augmented by Caputo-type fractional memory integrals. The governing equations retain the classical second-order time derivative as the leading operator and incorporate weakly singular convolution integrals with power-law kernels t−α, α∈(0,1), modelling viscoelastic memory damping rather than replacing the wave operator. The spatial discretisation employs Chebyshev–Gauss–Lobatto collocation, while the temporal integration uses a Newmark scheme (βNM=1/4) with the spatial operator treated implicitly and both the L1 memory sums and the cubic nonlinearities evaluated explicitly at the known time level; a linear extrapolation of the nonlinear terms eliminates the need for Newton–Raphson iterations. The disparate memory tails arising from two distinct fractional orders α≠β are compressed by independent Sum-of-Exponentials (SOE) approximations, reducing the per-step memory cost from O(Nt) to O(p+Nexp) and the total complexity from O(Nt2) to O(Nt(p+Nexp)). A rigorous stability estimate and a global convergence bound are established using a discrete Gronwall inequality. Numerical experiments confirm the temporal convergence rate O(Δtmin(2−α,2−β)), spectral spatial accuracy, and the practical speedup afforded by the SOE acceleration. A solitary wave collision scenario illustrates the method’s capability to capture asymmetric dispersive wakes generated by the fractional memory. The algorithmic architecture is dimension-independent by construction; a concrete extension pathway to multi-dimensional tensor-product Chebyshev grids, including Kronecker-product operators and Sylvester-based solvers, is presented.

Abstract: Numerical simulation of the multi-dimensional space-fractional Cahn-Hilliard equation faces two main computational challenges: the inherent temporal accuracy limitations of standard scalar auxiliary variable (SAV) methods and the escalating computational cost in high-dimensional domains. To address these issues, this study constructs a fully discrete algorithmic framework integrating a second-order backward differentiation formula (SAV-BDF2) with a sixth-order centered difference scheme. Under this formulation, we rigorously prove unconditional energy stability and establish the theoretical validity of the dual temporal and spatial accuracy. To solve the resulting indefinite algebraic systems, a minimal residual solver is paired with a sine-transform block diagonal preconditioner. Additionally, a hardware-level Vectorized Tensor Processing (VTP) architecture is deployed to resolve cache thrashing caused by non-contiguous memory access during multidimensional tensor evaluations. Numerical experiments in 3D to 8D domains demonstrate that the framework improves memory throughput and reduces execution time. By avoiding standard hardware execution inefficiencies, this integrated strategy provides an efficient numerical solution for large-scale simulations of high-dimensional fractional systems.