qualia_core_db/solvers/special_functions/
airy.rs

1//! Airy functions `Ai(x)` and `Bi(x)` via their Maclaurin series. Accurate for moderate
2//! `|x|` (the series converge for all `x` but lose digits for large argument).
3//!
4//! `Ai = α·f − β·g`, `Bi = √3·(α·f + β·g)`, where `α = Ai(0) = 3^{-2/3}/Γ(2/3)`,
5//! `β = −Ai'(0) = 3^{-1/3}/Γ(1/3)`, and `f`, `g` are the two hypergeometric series.
6
7const ALPHA: f64 = 0.355_028_053_887_817_24; // Ai(0)
8const BETA: f64 = 0.258_819_403_792_806_8; // −Ai'(0)
9const MAX_TERMS: usize = 200;
10
11/// `f(x) = Σ_k [∏(3j+1)] x^{3k}/(3k)!`, ratio `t_k/t_{k-1} = x³/((3k)(3k−1))`.
12fn series_f(x: f64) -> f64 {
13    let x3 = x * x * x;
14    let mut t = 1.0;
15    let mut sum = 0.0;
16    for k in 0..MAX_TERMS {
17        sum += t;
18        let k1 = (k + 1) as f64;
19        t *= x3 / ((3.0 * k1) * (3.0 * k1 - 1.0));
20        if t.abs() < 1e-18 {
21            break;
22        }
23    }
24    sum
25}
26
27/// `g(x) = Σ_k [∏(3j+2)] x^{3k+1}/(3k+1)!`, ratio `t_k/t_{k-1} = x³/((3k+1)(3k))`.
28fn series_g(x: f64) -> f64 {
29    let x3 = x * x * x;
30    let mut t = x;
31    let mut sum = 0.0;
32    for k in 0..MAX_TERMS {
33        sum += t;
34        let k1 = (k + 1) as f64;
35        t *= x3 / ((3.0 * k1 + 1.0) * (3.0 * k1));
36        if t.abs() < 1e-18 {
37            break;
38        }
39    }
40    sum
41}
42
43/// Airy function of the first kind `Ai(x)`.
44pub fn airy_ai(x: f64) -> f64 {
45    ALPHA * series_f(x) - BETA * series_g(x)
46}
47
48/// Airy function of the second kind `Bi(x)`.
49pub fn airy_bi(x: f64) -> f64 {
50    3.0_f64.sqrt() * (ALPHA * series_f(x) + BETA * series_g(x))
51}
52
53#[cfg(test)]
54mod tests {
55    use super::*;
56    const TOL: f64 = 1e-9;
57
58    #[test]
59    fn airy_at_origin() {
60        assert!((airy_ai(0.0) - ALPHA).abs() < TOL);
61        assert!((airy_bi(0.0) - 3.0_f64.sqrt() * ALPHA).abs() < TOL);
62    }
63
64    #[test]
65    fn airy_table_values() {
66        assert!((airy_ai(1.0) - 0.135_292_416_312_881_4).abs() < 1e-7);
67        assert!((airy_bi(1.0) - 1.207_423_594_952_871_3).abs() < 1e-7);
68        assert!((airy_ai(-1.0) - 0.535_560_883_292_352_2).abs() < 1e-7);
69    }
70}
qualia_core_db/solvers/special_functions/airy.rs

qualia_core_db/solvers/special_functions/
airy.rs
