We propose an efficient high-order quadrature scheme for approximating convolution integrals with weakly singular kernels, specifically \( g_\alpha (|x-y|) = |x-y|^{-\alpha } \ \text {for } 0< \alpha < 1, \ g_\alpha (|x-y|) = \log |x-y| \quad \text {for } \alpha = 0. \) The method partitions the domain into subdomains (or patches) and analytically resolves the singularity using a polynomial change of variables centered at the target point. A central theoretical contribution of this work is the precise characterization of convergence rates in terms of the kernel parameter \(\alpha \) and the transformation degree p, establishing that high-order accuracy is achieved whenever \(p(1-\alpha ) \in \mathbb {N}\) . Moreover, we prove that splitting the integral at the singularity leads to a doubling of the convergence order compared to the unsplit case. The scheme extends naturally to bounded curves in two dimensions for arbitrary weakly singular kernels. For computational efficiency, the domain is refined into multiple subdomains while maintaining a fixed number of quadrature nodes per subdomain, enabling the use of fast summation techniques to accelerate the evaluation of regular integrals. Numerical experiments validate the theoretical findings, consistently achieving machine precision. Finally, we demonstrate the applicability of the method by solving two-dimensional surface wave scattering problems on complex geometries, attaining dispersion-free, high-accuracy solutions that scale efficiently to large problem sizes.
