Saddle-point adapted quadrature for canonical partition function integrals #62
Description
Problem
The canonical partition function integrals in ThermalModelCanonical::CalculatePartitionFunctions() become expensive at large volumes because the integrand develops a sharp peak around φ = 0 whose width scales as ~1/√V, requiring more quadrature points to resolve.
Idea
Factor out the peak analytically. The integrand is:
f(φ) = exp(Σᵢ Nᵢ [cos(qᵢ·φ) - 1]) · cos(Q·φ)
Expanding cos(qᵢ·φ) - 1 ≈ -½(qᵢ·φ)² near the peak gives a Gaussian approximation with curvature matrix A = Σᵢ Nᵢ qᵢ qᵢᵀ, computable from the GCE multiplicities Nsx[i] already available. Rewrite as:
f(φ) = exp(-½ φᵀ A φ) · g(φ)
where g(φ) captures only deviations from Gaussianity. Handle the Gaussian weight analytically; apply the trapezoidal rule to the flat remainder g(φ).
Properties
- No equation solving needed — saddle point is always at φ = 0, curvature matrix comes from existing GCE multiplicities
- Convergence improves with volume — at large V, g(φ) → cos(Q·φ), so fewer points suffice
- Presumably graceful small-V limit — A → 0, Gaussian → 1, reduces to the plain trapezoidal rule
Optional refinement
For large non-zero Q, we will still get large oscillatory behavior from cos(Q·φ).
To address that, shift the expansion from φ = 0 to the true saddle point φ* of Σᵢ Nᵢ [cos(qᵢ·φ) - 1] - Q·φ, making g(φ) even flatter for large Q. Finding φ* is equivalent to solving for GCE chemical potentials.