Eigenvalue-stabilized high-order Gauss-Laguerre quadrature with practical a posteriori error estimation for double integrals in the average kernel method

The accurate evaluation of double integrals with exponential-decay kernels is central to the average kernel (AK) method for solving the Smoluchowski coagulation equation (SCE), a cornerstone model in aerosol physics, combustion, and astrophysical dust evolution. However, its precision has been fundamentally limited by two intertwined challenges: the numerical instability of high-order Gauss-Laguerre quadrature and the lack of a practical error control mechanism for complex kernels. This paper introduces a dual-faceted numerical framework that addresses both issues, thereby unlocking the full predictive potential of the AK method for high-fidelity simulations of coagulation-driven systems. First, we reformulate the computation of quadrature nodes and weights as a symmetric eigenvalue problem, enabling stable and reliable calculations at orders up to n > 10 4 , bypassing the severe ill-conditioning of traditional root-finding methods. Second, and more importantly, we develop a novel a posteriori error estimation technique that does not require derivatives of the kernel function. By analyzing the decay of the iterative error ϵ n , n = | Q n + 1 , n + 1 − Q n , n | , we derive a robust power-law model ( ϵ n , n ∝ n C ) that accurately predicts convergence rates and truncation errors. This model serves as a practical surrogate for classical derivative-based bounds and provides a computationally cheap mechanism for error control, a critical need in time-dependent multiphysics simulations. Numerical experiments demonstrate that our framework achieves a reduction in relative error by 2–3 orders of magnitude (factors of 150 to 250) while preserving the homogeneity of physical kernels, thereby resolving long-standing paradoxes like the zero-sedimentation kernel and ensuring physically consistent predictions. The proposed methodology is not only transformative for coagulation simulations but also serves as a general-purpose tool for high-fidelity integration over semi-infinite domains with exponential weights, with direct applications in radiation transport, stochastic process integration, and uncertainty quantification.