Non‑Polynomial Numerical Quadrature via Higher‑Order Osculating Conics with Operator‑Based Jet Extraction from Discrete Invariant Projection Spaces

We develop a composite numerical quadrature scheme in which the integrand is locally approximated by arcs of osculating conics—circles or ellipses—rather than by polynomial interpolants. In each subinterval the arc matches the function value, slope, and curvature (and optionally higher derivatives) at the midpoint. The area under the arc is computed analytically via a closed form that subtracts the conic segment from its bounding rectangle and adds a rectangular support. We provide a complete error analysis: for the composite circle method (MCOC) we prove global convergence of order O(h⁴) under the hypothesis f ∈ C⁶, and we give the exact leading error constant. A sharp stability condition determines exactly when the arc remains real on the subinterval. The extension to fourth‑order osculating ellipses (MEO4) yields O(h⁶) convergence; a classification theorem guarantees fallback to MCOC when the conic is not an ellipse. For integration over a torus, we prove that mixed derivative terms vanish on a symmetric Cartesian mesh, preserving separability without grid rotation. To eliminate the need for costly analytical derivatives, we link the method to Discrete Invariant Projection Spaces (DIPS). Exploiting the Bose–Mesner algebra of the sampling mesh, we construct a global deconvolution operator that recovers the local differential jet from a single function evaluation per node, and we prove that the quadrature order is preserved. We give a full complexity and conditioning analysis, and we establish consistency, stability, and convergence theorems in the sense of Lax–Richtmyer. The resulting scheme requires the same number of function evaluations as a Riemann sum while delivering high order and geometric adaptivity. All proofs are complete and explicit; the critical algebraic calculations—the spectral deconvolution kernel, the error splitting in jet extraction, the exact cancellation of torsion terms on the torus, and the sixth derivative of the osculating ellipse—are presented in full detail. Numerical experiments on transcendental and algebraic test functions confirm the theoretical orders and demonstrate that, for functions with smooth curvature, MCOC can surpass Simpson’s rule with the same number of subintervals.