1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
//! Braess-Hackbusch exponential sum approximation of 1/r.
//!
//! The Coulomb/Newton kernel 1/r is approximated as a sum of Gaussians:
//!
//! ```text
//! 1/r ≈ (2/√π) Σ_{k=1}^{R_G} c_k exp(-α_k r²)
//! ```
//!
//! via discretization of the Gaussian integral representation
//! `1/r = (2/√π) ∫₀^∞ exp(-t² r²) dt` using the substitution `t = exp(u)`
//! and the trapezoidal rule.
//!
//! Reference: Braess & Hackbusch, "Approximation of 1/x by exponential sums
//! in [1,∞)", IMA J. Numer. Anal. 25(4), 685–697, 2005.
/// Exponential sum coefficients for 1/r ≈ Σ c_k exp(-α_k r²).
pub struct ExponentialSumCoefficients {
/// Weights c_k (include the 2/√π prefactor and step h).
pub c: Vec<f64>,
/// Exponents α_k = exp(2 u_k).
pub alpha: Vec<f64>,
/// Number of terms R_G.
pub r_g: usize,
}
impl ExponentialSumCoefficients {
/// Compute exponential sum coefficients for approximating 1/r
/// on [delta, r_max] to accuracy epsilon.
///
/// Uses the Braess-Hackbusch formula: substitute t = exp(u) in the
/// Gaussian integral 1/r = (2/√π) ∫ exp(-t²r²) dt, then apply
/// the trapezoidal rule with optimized step h.
pub fn compute(delta: f64, r_max: f64, epsilon: f64) -> Self {
assert!(delta > 0.0, "delta must be positive");
assert!(r_max > delta, "r_max must exceed delta");
assert!(epsilon > 0.0 && epsilon < 1.0, "epsilon must be in (0,1)");
// The BH formula gives R_G = O(log(1/ε) · log(R/δ)) terms.
// Optimal step h ≈ π / (2 · √(log(1/ε) · log(R/δ))) balances the
// trapezoidal discretization error against the support width.
let ln_r = (r_max / delta).ln().max(1.0);
let ln_eps = (1.0 / epsilon).ln().max(1.0);
let h = std::f64::consts::PI / (2.0 * (ln_eps * ln_r).sqrt());
// Cutoffs in u-space (unnormalized, i.e., for r in [delta, r_max]):
// Upper: exp(-exp(2u_max) * delta²) < machine_eps
// → exp(2u_max) > 36 / delta² → u_max = 0.5 * ln(36 / delta²)
let u_max = 0.5 * (36.0 / (delta * delta)).ln();
// Lower: for u → -∞, exp(2u) → 0 so the Gaussian ≈ 1, and the
// integrand ≈ exp(u). The tail integral from -∞ to u_min is ≈ exp(u_min).
// The full integral is √π/(2r). So relative error at r_max is
// 2·r_max·exp(u_min)/√π. For this < ε:
// u_min < ln(ε·√π / (2·r_max))
let u_min = (epsilon * std::f64::consts::PI.sqrt() / (2.0 * r_max)).ln();
let two_over_sqrt_pi = std::f64::consts::FRAC_2_SQRT_PI;
let n_neg = ((0.0 - u_min) / h).ceil() as isize;
let n_pos = (u_max / h).ceil() as isize;
let mut c = Vec::with_capacity((n_neg + n_pos + 1) as usize);
let mut alpha = Vec::with_capacity((n_neg + n_pos + 1) as usize);
for k in -n_neg..=n_pos {
let u = k as f64 * h;
let exp_u = u.exp();
let alpha_k = exp_u * exp_u; // exp(2u)
let c_k = two_over_sqrt_pi * h * exp_u;
// Skip negligible weights
if c_k < 1e-20 {
continue;
}
// Skip terms where the Gaussian is dead at r = delta
if (-alpha_k * delta * delta).exp() < 1e-18 {
continue;
}
c.push(c_k);
alpha.push(alpha_k);
}
let r_g = c.len();
Self { c, alpha, r_g }
}
/// Evaluate the exponential sum approximation at distance r.
pub fn evaluate(&self, r: f64) -> f64 {
let r2 = r * r;
let mut sum = 0.0;
for k in 0..self.r_g {
sum += self.c[k] * (-self.alpha[k] * r2).exp();
}
sum
}
/// Absolute error |approx(r) - 1/r|.
pub fn error_at(&self, r: f64) -> f64 {
(self.evaluate(r) - 1.0 / r).abs()
}
/// Relative error |approx(r)*r - 1|.
pub fn relative_error_at(&self, r: f64) -> f64 {
(self.evaluate(r) * r - 1.0).abs()
}
/// Near-field correction: the difference between the exact 1/r and the
/// Gaussian sum approximation, integrated over a small ball of radius delta.
///
/// Returns the correction potential at a point where the density is `rho`:
/// Φ_corr = -G * rho * ∫_{B_δ} [1/|y| - GS(|y|)] dy
///
/// where GS is the Gaussian sum approximation.
/// The integral is computed analytically using the difference:
/// ∫_0^δ [1/r - GS(r)] * 4πr² dr
///
/// This 0th-order correction (Exl, Mauser & Zhang, JCP 2016) handles
/// the singularity mismatch at r < δ.
pub fn near_field_correction_integral(&self, delta: f64) -> f64 {
let n_quad = 200;
let dr = delta / n_quad as f64;
let mut integral = 0.0;
for i in 0..n_quad {
let r = (i as f64 + 0.5) * dr;
let exact = 1.0 / r;
let approx = self.evaluate(r);
let diff = exact - approx;
integral += diff * 4.0 * std::f64::consts::PI * r * r * dr;
}
integral
}
/// Second-order near-field correction integrals (Exl, Mauser & Zhang 2016).
///
/// Returns `(I_0, I_2)` where:
/// - `I_0 = ∫₀^δ [1/r - GS(r)] 4πr² dr` (0th-order, same as `near_field_correction_integral`)
/// - `I_2 = ∫₀^δ [1/r - GS(r)] 4πr⁴/6 dr` (2nd-order, couples to ∇²ρ)
///
/// The second-order correction is: `Φ_corr_2(x) = -G · I_2 · ∇²ρ(x)`
pub fn near_field_correction_second_order(&self, delta: f64) -> (f64, f64) {
let n_quad = 200;
let dr = delta / n_quad as f64;
let mut i0 = 0.0;
let mut i2 = 0.0;
for i in 0..n_quad {
let r = (i as f64 + 0.5) * dr;
let exact = 1.0 / r;
let approx = self.evaluate(r);
let diff = exact - approx;
let r2 = r * r;
i0 += diff * 4.0 * std::f64::consts::PI * r2 * dr;
i2 += diff * 4.0 * std::f64::consts::PI * r2 * r2 / 6.0 * dr;
}
(i0, i2)
}
/// Maximum relative error on [delta, r_max] sampled at n_test log-spaced points.
pub fn max_relative_error(&self, delta: f64, r_max: f64, n_test: usize) -> f64 {
let mut max_err = 0.0f64;
for i in 0..n_test {
let t = i as f64 / (n_test - 1) as f64;
let r = delta * (r_max / delta).powf(t);
max_err = max_err.max(self.relative_error_at(r));
}
max_err
}
}
#[cfg(test)]
mod tests {
use super::*;
#[test]
fn exp_sum_accuracy_1e4() {
let delta = 0.1;
let r_max = 10.0;
let epsilon = 1e-4;
let coeffs = ExponentialSumCoefficients::compute(delta, r_max, epsilon);
// R_G should scale as O(log(1/ε) * log(R/δ))
// For ε=1e-4, R/δ=100: expect ~30-150 terms
assert!(coeffs.r_g > 5, "Too few terms: {}", coeffs.r_g);
let max_err = coeffs.max_relative_error(delta, r_max, 500);
assert!(
max_err < epsilon,
"Max relative error {max_err} exceeds tolerance {epsilon}, R_G={}",
coeffs.r_g
);
}
#[test]
fn exp_sum_accuracy_1e8() {
let delta = 0.05;
let r_max = 20.0;
let epsilon = 1e-8;
let coeffs = ExponentialSumCoefficients::compute(delta, r_max, epsilon);
let max_err = coeffs.max_relative_error(delta, r_max, 500);
assert!(
max_err < epsilon,
"Max relative error {max_err} exceeds tolerance {epsilon}, R_G={}",
coeffs.r_g
);
}
#[test]
fn exp_sum_evaluate_at_one() {
let coeffs = ExponentialSumCoefficients::compute(0.1, 10.0, 1e-6);
let val = coeffs.evaluate(1.0);
assert!(
(val - 1.0).abs() < 1e-5,
"evaluate(1.0) = {val}, expected ~1.0"
);
}
#[test]
fn exp_sum_scaling() {
// Smaller step h (tighter tolerance) should give more terms
let c1 = ExponentialSumCoefficients::compute(0.1, 10.0, 1e-4);
let c2 = ExponentialSumCoefficients::compute(0.1, 10.0, 1e-8);
assert!(
c2.r_g > c1.r_g,
"Higher accuracy should require more terms: {} vs {}",
c2.r_g,
c1.r_g
);
}
}