qualia_core_db/solvers/special_functions/
bessel.rs

1//! Bessel functions of integer order: `J_n`, `Y_n` (first/second kind) and the modified
2//! `I_n`, `K_n`. The order-0 functions come from their convergent power series (with the
3//! log + harmonic-number terms for the second kinds); order 1 from the Wronskian
4//! relations; higher orders from the standard upward recurrences. Accurate for moderate
5//! `|x|` (the series converge for all `x` but lose digits for large argument — documented).
6
7const EULER_GAMMA: f64 = 0.577_215_664_901_532_9;
8const MAX_TERMS: usize = 400;
9
10fn factorial_f(n: u32) -> f64 {
11    (1..=n).fold(1.0, |a, k| a * k as f64)
12}
13
14/// `H_k = 1 + 1/2 + … + 1/k`.
15fn harmonic(k: u32) -> f64 {
16    (1..=k).fold(0.0, |a, j| a + 1.0 / j as f64)
17}
18
19fn j_nonneg(n: u32, x: f64) -> f64 {
20    // Σ_{m≥0} (−1)^m (x/2)^{2m+n} / (m! (m+n)!)
21    let h2 = (x / 2.0) * (x / 2.0);
22    let mut a = (x / 2.0).powi(n as i32) / factorial_f(n); // m = 0
23    let mut sum = 0.0;
24    let mut sign = 1.0;
25    for m in 0..MAX_TERMS {
26        sum += sign * a;
27        a *= h2 / (((m + 1) as f64) * ((m + 1 + n as usize) as f64));
28        sign = -sign;
29        if a.abs() < 1e-18 && m as u32 > n {
30            break;
31        }
32    }
33    sum
34}
35
36fn i_nonneg(n: u32, x: f64) -> f64 {
37    // Σ_{m≥0} (x/2)^{2m+n} / (m! (m+n)!)  — same as J but all-positive.
38    let h2 = (x / 2.0) * (x / 2.0);
39    let mut a = (x / 2.0).powi(n as i32) / factorial_f(n);
40    let mut sum = 0.0;
41    for m in 0..MAX_TERMS {
42        sum += a;
43        a *= h2 / (((m + 1) as f64) * ((m + 1 + n as usize) as f64));
44        if a.abs() < 1e-18 && m as u32 > n {
45            break;
46        }
47    }
48    sum
49}
50
51/// Bessel function of the first kind `J_n(x)`, integer order (any sign). Defined for all
52/// real `x`. `J_{-n} = (−1)^n J_n`.
53pub fn bessel_j(n: i32, x: f64) -> f64 {
54    let m = n.unsigned_abs();
55    let v = j_nonneg(m, x);
56    if n < 0 && m % 2 == 1 {
57        -v
58    } else {
59        v
60    }
61}
62
63/// Modified Bessel function of the first kind `I_n(x)`, integer order. `I_{-n} = I_n`.
64pub fn bessel_i(n: i32, x: f64) -> f64 {
65    i_nonneg(n.unsigned_abs(), x)
66}
67
68fn y0(x: f64) -> f64 {
69    // Y_0 = (2/π)(ln(x/2)+γ)J_0 + (2/π) Σ_{k≥1} (−1)^{k+1} H_k/(k!)^2 (x/2)^{2k}
70    let h2 = (x / 2.0) * (x / 2.0);
71    let mut series = 0.0;
72    let mut term = 1.0; // (x/2)^{2k}/(k!)^2 accumulator, starts at k=0 value 1
73    let mut sign = 1.0; // (−1)^{k+1} for k=1 is +1
74    for k in 1..MAX_TERMS as u32 {
75        term *= h2 / (k as f64 * k as f64); // now (x/2)^{2k}/(k!)^2
76        series += sign * harmonic(k) * term;
77        sign = -sign;
78        if term * harmonic(k) < 1e-18 {
79            break;
80        }
81    }
82    let two_over_pi = 2.0 / core::f64::consts::PI;
83    two_over_pi * ((x / 2.0).ln() + EULER_GAMMA) * j_nonneg(0, x) + two_over_pi * series
84}
85
86fn k0(x: f64) -> f64 {
87    // K_0 = −(ln(x/2)+γ)I_0 + Σ_{k≥1} H_k/(k!)^2 (x/2)^{2k}
88    let h2 = (x / 2.0) * (x / 2.0);
89    let mut series = 0.0;
90    let mut term = 1.0;
91    for k in 1..MAX_TERMS as u32 {
92        term *= h2 / (k as f64 * k as f64);
93        series += harmonic(k) * term;
94        if term * harmonic(k) < 1e-18 {
95            break;
96        }
97    }
98    -((x / 2.0).ln() + EULER_GAMMA) * i_nonneg(0, x) + series
99}
100
101/// Bessel function of the second kind `Y_n(x)`, integer order `n ≥ 0`. Requires `x > 0`
102/// (singular at the origin) → `None` otherwise. Order 1 via the Wronskian
103/// `J_1 Y_0 − J_0 Y_1 = 2/(πx)`; higher via `Y_{n+1} = (2n/x)Y_n − Y_{n-1}`.
104pub fn bessel_y(n: u32, x: f64) -> Option<f64> {
105    if x <= 0.0 {
106        return None;
107    }
108    let y0v = y0(x);
109    if n == 0 {
110        return Some(y0v);
111    }
112    let j0 = j_nonneg(0, x);
113    if j0.abs() < 1e-300 {
114        return None; // at a zero of J_0 the Wronskian solve is ill-posed
115    }
116    let y1 = (bessel_j(1, x) * y0v - 2.0 / (core::f64::consts::PI * x)) / j0;
117    if n == 1 {
118        return Some(y1);
119    }
120    let (mut ym1, mut yn) = (y0v, y1);
121    for k in 1..n {
122        let ynext = (2.0 * k as f64 / x) * yn - ym1;
123        ym1 = yn;
124        yn = ynext;
125    }
126    Some(yn)
127}
128
129/// Modified Bessel function of the second kind `K_n(x)`, integer order `n ≥ 0`. Requires
130/// `x > 0` → `None` otherwise. Order 1 via the Wronskian `I_0 K_1 + I_1 K_0 = 1/x`;
131/// higher via `K_{n+1} = (2n/x)K_n + K_{n-1}`.
132pub fn bessel_k(n: u32, x: f64) -> Option<f64> {
133    if x <= 0.0 {
134        return None;
135    }
136    let k0v = k0(x);
137    if n == 0 {
138        return Some(k0v);
139    }
140    let i0 = i_nonneg(0, x);
141    let k1 = (1.0 / x - bessel_i(1, x) * k0v) / i0;
142    if n == 1 {
143        return Some(k1);
144    }
145    let (mut km1, mut kn) = (k0v, k1);
146    for k in 1..n {
147        let knext = (2.0 * k as f64 / x) * kn + km1;
148        km1 = kn;
149        kn = knext;
150    }
151    Some(kn)
152}
153
154#[cfg(test)]
155mod tests {
156    use super::*;
157    const TOL: f64 = 1e-7;
158
159    #[test]
160    fn first_kind_table() {
161        assert!((bessel_j(0, 0.0) - 1.0).abs() < TOL);
162        assert!((bessel_j(1, 0.0)).abs() < TOL);
163        assert!((bessel_j(0, 1.0) - 0.765_197_686_557_966_5).abs() < TOL);
164        assert!((bessel_j(1, 1.0) - 0.440_050_585_744_933_5).abs() < TOL);
165        assert!((bessel_j(2, 2.0) - 0.352_834_028_615_815_5).abs() < TOL);
166        // J_{-1} = −J_1
167        assert!((bessel_j(-1, 1.0) + bessel_j(1, 1.0)).abs() < TOL);
168    }
169
170    #[test]
171    fn modified_first_kind_table() {
172        assert!((bessel_i(0, 0.0) - 1.0).abs() < TOL);
173        assert!((bessel_i(0, 1.0) - 1.266_065_877_752_008_4).abs() < TOL);
174        assert!((bessel_i(1, 1.0) - 0.565_159_103_992_485_0).abs() < TOL);
175    }
176
177    #[test]
178    fn second_kind_table_and_domain() {
179        assert!(bessel_y(0, -1.0).is_none()); // x ≤ 0 fails closed
180        assert!((bessel_y(0, 1.0).unwrap() - 0.088_256_964_215_676_96).abs() < TOL);
181        assert!((bessel_y(1, 1.0).unwrap() + 0.781_212_821_300_288_7).abs() < TOL);
182        assert!((bessel_k(0, 1.0).unwrap() - 0.421_024_438_240_708_3).abs() < TOL);
183        assert!((bessel_k(1, 1.0).unwrap() - 0.601_907_230_197_234_6).abs() < TOL);
184        assert!(bessel_k(2, 0.0).is_none());
185    }
186}
qualia_core_db/solvers/special_functions/bessel.rs

qualia_core_db/solvers/special_functions/
bessel.rs
