1pub fn adaptive_simpson(f: impl Fn(f64) -> f64, a: f64, b: f64, tol: f64, max_depth: usize) -> f64 {
2 adaptive_simpson_rec(&f, a, b, tol, max_depth, 0)
3}
4
5fn adaptive_simpson_rec(
6 f: &impl Fn(f64) -> f64,
7 a: f64,
8 b: f64,
9 tol: f64,
10 max_depth: usize,
11 depth: usize,
12) -> f64 {
13 let mid = 0.5 * (a + b);
14 let h = b - a;
15 let s_whole = h / 6.0 * (f(a) + 4.0 * f(mid) + f(b));
16 let s_left = (mid - a) / 6.0 * (f(a) + 4.0 * f(0.5 * (a + mid)) + f(mid));
17 let s_right = (b - mid) / 6.0 * (f(mid) + 4.0 * f(0.5 * (mid + b)) + f(b));
18 let s_split = s_left + s_right;
19
20 if depth >= max_depth || (s_split - s_whole).abs() < 15.0 * tol {
21 return s_split + (s_split - s_whole) / 15.0;
22 }
23
24 adaptive_simpson_rec(f, a, mid, tol * 0.5, max_depth, depth + 1)
25 + adaptive_simpson_rec(f, mid, b, tol * 0.5, max_depth, depth + 1)
26}
27
28pub fn adaptive_trapezoid(
29 f: impl Fn(f64) -> f64,
30 a: f64,
31 b: f64,
32 tol: f64,
33 max_depth: usize,
34) -> f64 {
35 adaptive_trapezoid_rec(&f, a, b, tol, max_depth, 0)
36}
37
38fn adaptive_trapezoid_rec(
39 f: &impl Fn(f64) -> f64,
40 a: f64,
41 b: f64,
42 tol: f64,
43 max_depth: usize,
44 depth: usize,
45) -> f64 {
46 let mid = 0.5 * (a + b);
47 let h = b - a;
48 let whole = 0.5 * h * (f(a) + f(b));
49 let left = 0.5 * (mid - a) * (f(a) + f(mid));
50 let right = 0.5 * (b - mid) * (f(mid) + f(b));
51 let split = left + right;
52
53 if depth >= max_depth || (split - whole).abs() < 3.0 * tol {
54 return split + (split - whole) / 3.0;
55 }
56
57 adaptive_trapezoid_rec(f, a, mid, tol * 0.5, max_depth, depth + 1)
58 + adaptive_trapezoid_rec(f, mid, b, tol * 0.5, max_depth, depth + 1)
59}
60
61pub fn improper_integral_transform(f: impl Fn(f64) -> f64, a: f64, n: usize) -> f64 {
62 let g = |t: f64| {
63 if t.abs() < 1e-30 {
64 return 0.0;
65 }
66 let x = a + (1.0 - t) / t;
67 f(x) / (t * t)
68 };
69 super::quadrature::simpson(g, 1e-10, 1.0, n)
70}
71
72pub fn adaptive_gauss_kronrod(
73 f: impl Fn(f64) -> f64,
74 a: f64,
75 b: f64,
76 tol: f64,
77 max_depth: usize,
78) -> f64 {
79 adaptive_gk_rec(&f, a, b, tol, max_depth, 0)
80}
81
82fn adaptive_gk_rec(
83 f: &impl Fn(f64) -> f64,
84 a: f64,
85 b: f64,
86 tol: f64,
87 max_depth: usize,
88 depth: usize,
89) -> f64 {
90 let (val, err) = super::quadrature::gauss_kronrod_15(f, a, b);
91 if depth >= max_depth || err < tol {
92 return val;
93 }
94 let mid = 0.5 * (a + b);
95 adaptive_gk_rec(f, a, mid, tol * 0.5, max_depth, depth + 1)
96 + adaptive_gk_rec(f, mid, b, tol * 0.5, max_depth, depth + 1)
97}
98
99pub fn adaptive_midpoint(
100 f: impl Fn(f64) -> f64,
101 a: f64,
102 b: f64,
103 tol: f64,
104 max_depth: usize,
105) -> f64 {
106 adaptive_midpoint_rec(&f, a, b, tol, max_depth, 0)
107}
108
109fn adaptive_midpoint_rec(
110 f: &impl Fn(f64) -> f64,
111 a: f64,
112 b: f64,
113 tol: f64,
114 max_depth: usize,
115 depth: usize,
116) -> f64 {
117 let mid = 0.5 * (a + b);
118 let h = b - a;
119 let whole = h * f(mid);
120 let left = 0.5 * h * f(0.5 * (a + mid));
121 let right = 0.5 * h * f(0.5 * (mid + b));
122 let split = left + right;
123 if depth >= max_depth || (split - whole).abs() < 3.0 * tol {
124 return split;
125 }
126 adaptive_midpoint_rec(f, a, mid, tol * 0.5, max_depth, depth + 1)
127 + adaptive_midpoint_rec(f, mid, b, tol * 0.5, max_depth, depth + 1)
128}
129
130pub fn adaptive_boole(f: impl Fn(f64) -> f64, a: f64, b: f64, tol: f64, max_depth: usize) -> f64 {
131 adaptive_boole_rec(&f, a, b, tol, max_depth, 0)
132}
133
134fn adaptive_boole_rec(
135 f: &impl Fn(f64) -> f64,
136 a: f64,
137 b: f64,
138 tol: f64,
139 max_depth: usize,
140 depth: usize,
141) -> f64 {
142 let h = b - a;
143 let x0 = a;
144 let x1 = a + 0.25 * h;
145 let x2 = a + 0.5 * h;
146 let x3 = a + 0.75 * h;
147 let x4 = b;
148 let whole = h / 90.0 * (7.0 * f(x0) + 32.0 * f(x1) + 12.0 * f(x2) + 32.0 * f(x3) + 7.0 * f(x4));
149 let mid = 0.5 * (a + b);
150 let left = super::quadrature::simpson(f, a, mid, 4);
151 let right = super::quadrature::simpson(f, mid, b, 4);
152 let split = left + right;
153 if depth >= max_depth || (split - whole).abs() < 15.0 * tol {
154 return split + (split - whole) / 15.0;
155 }
156 adaptive_boole_rec(f, a, mid, tol * 0.5, max_depth, depth + 1)
157 + adaptive_boole_rec(f, mid, b, tol * 0.5, max_depth, depth + 1)
158}
159
160pub fn double_exponential(f: impl Fn(f64) -> f64, a: f64, b: f64, n: usize) -> f64 {
161 let pi_half = std::f64::consts::FRAC_PI_2;
162 let mid = 0.5 * (a + b);
163 let half = 0.5 * (b - a);
164 let h = 6.0 / n as f64;
165 let mut sum = 0.0;
166 for k in 0..=n {
167 let t = -3.0 + k as f64 * h;
168 let phi = pi_half * t.sinh();
169 let phi_prime = pi_half * t.cosh();
170 let x_t = phi.tanh();
171 let cosh_phi = phi.cosh();
172 let w = phi_prime / (cosh_phi * cosh_phi);
173 if w.abs() < 1e-30 {
174 continue;
175 }
176 let xm = mid + half * x_t;
177 if xm > a && xm < b {
178 sum += w * f(xm);
179 }
180 }
181 sum * half * h
182}
183
184pub fn cauchy_principal_value(
185 f: impl Fn(f64) -> f64,
186 a: f64,
187 b: f64,
188 singularity: f64,
189 epsilon: f64,
190 n: usize,
191) -> f64 {
192 let left = if singularity - epsilon > a {
193 super::quadrature::simpson(&f, a, singularity - epsilon, n)
194 } else {
195 0.0
196 };
197 let right = if singularity + epsilon < b {
198 super::quadrature::simpson(&f, singularity + epsilon, b, n)
199 } else {
200 0.0
201 };
202 left + right
203}
204
205pub fn improper_both_infinite(f: impl Fn(f64) -> f64, n: usize) -> f64 {
206 let g = |t: f64| {
207 let x = t / (1.0 - t * t);
208 let dx = (1.0 + t * t) / ((1.0 - t * t) * (1.0 - t * t));
209 f(x) * dx
210 };
211 super::quadrature::simpson(g, -1.0 + 1e-10, 1.0 - 1e-10, n)
212}
213
214pub fn improper_left_infinite(f: impl Fn(f64) -> f64, b: f64, n: usize) -> f64 {
215 let g = |t: f64| {
216 if t.abs() < 1e-30 {
217 return 0.0;
218 }
219 let x = b - (1.0 - t) / t;
220 f(x) / (t * t)
221 };
222 super::quadrature::simpson(g, 1e-10, 1.0, n)
223}
224
225pub fn numerical_derivative_central(f: impl Fn(f64) -> f64, x: f64, h: f64) -> f64 {
226 (f(x + h) - f(x - h)) / (2.0 * h)
227}
228
229pub fn numerical_second_derivative(f: impl Fn(f64) -> f64, x: f64, h: f64) -> f64 {
230 (f(x + h) - 2.0 * f(x) + f(x - h)) / (h * h)
231}
232
233pub fn numerical_derivative_5point(f: impl Fn(f64) -> f64, x: f64, h: f64) -> f64 {
234 (-f(x + 2.0 * h) + 8.0 * f(x + h) - 8.0 * f(x - h) + f(x - 2.0 * h)) / (12.0 * h)
235}
236
237pub fn numerical_integral_cumulative(
238 f: impl Fn(f64) -> f64,
239 a: f64,
240 b: f64,
241 n: usize,
242) -> Vec<(f64, f64)> {
243 let h = (b - a) / n as f64;
244 let mut result = Vec::with_capacity(n + 1);
245 result.push((a, 0.0));
246 let mut cumul = 0.0;
247 for i in 1..=n {
248 let x0 = a + (i - 1) as f64 * h;
249 let x1 = a + i as f64 * h;
250 cumul += 0.5 * h * (f(x0) + f(x1));
251 result.push((x1, cumul));
252 }
253 result
254}