proof_engine/solver/
spectral.rs

1//! Spectral methods — FFT-based PDE solving for periodic domains.
2
3/// 1D Discrete Fourier Transform (naive O(n²) — for correctness; production uses FFT).
4pub fn dft(input: &[(f64, f64)]) -> Vec<(f64, f64)> {
5    let n = input.len();
6    let mut output = vec![(0.0, 0.0); n];
7    for k in 0..n {
8        let mut re = 0.0;
9        let mut im = 0.0;
10        for j in 0..n {
11            let angle = -2.0 * std::f64::consts::PI * k as f64 * j as f64 / n as f64;
12            re += input[j].0 * angle.cos() - input[j].1 * angle.sin();
13            im += input[j].0 * angle.sin() + input[j].1 * angle.cos();
14        }
15        output[k] = (re, im);
16    }
17    output
18}
19
20/// 1D Inverse DFT.
21pub fn idft(input: &[(f64, f64)]) -> Vec<(f64, f64)> {
22    let n = input.len();
23    let mut output = vec![(0.0, 0.0); n];
24    for j in 0..n {
25        let mut re = 0.0;
26        let mut im = 0.0;
27        for k in 0..n {
28            let angle = 2.0 * std::f64::consts::PI * k as f64 * j as f64 / n as f64;
29            re += input[k].0 * angle.cos() - input[k].1 * angle.sin();
30            im += input[k].0 * angle.sin() + input[k].1 * angle.cos();
31        }
32        output[j] = (re / n as f64, im / n as f64);
33    }
34    output
35}
36
37/// Solve the 1D heat equation spectrally: u_t = α * u_xx on periodic domain [0, L].
38pub fn spectral_heat_step(u: &mut [f64], alpha: f64, dt: f64, dx: f64) {
39    let n = u.len();
40    let input: Vec<(f64, f64)> = u.iter().map(|&v| (v, 0.0)).collect();
41    let mut spectrum = dft(&input);
42
43    // Multiply by exp(-α * k² * dt) in frequency space
44    let l = n as f64 * dx;
45    for k in 0..n {
46        let kk = if k <= n / 2 { k as f64 } else { k as f64 - n as f64 };
47        let freq = 2.0 * std::f64::consts::PI * kk / l;
48        let decay = (-alpha * freq * freq * dt).exp();
49        spectrum[k].0 *= decay;
50        spectrum[k].1 *= decay;
51    }
52
53    let result = idft(&spectrum);
54    for i in 0..n { u[i] = result[i].0; }
55}
56
57/// Power spectrum (magnitude² of DFT coefficients).
58pub fn power_spectrum(signal: &[f64]) -> Vec<f64> {
59    let input: Vec<(f64, f64)> = signal.iter().map(|&v| (v, 0.0)).collect();
60    let spectrum = dft(&input);
61    spectrum.iter().map(|(re, im)| re * re + im * im).collect()
62}
63
64#[cfg(test)]
65mod tests {
66    use super::*;
67
68    #[test]
69    fn dft_idft_roundtrip() {
70        let input = vec![(1.0, 0.0), (2.0, 0.0), (3.0, 0.0), (4.0, 0.0)];
71        let spectrum = dft(&input);
72        let recovered = idft(&spectrum);
73        for (orig, rec) in input.iter().zip(recovered.iter()) {
74            assert!((orig.0 - rec.0).abs() < 1e-10, "orig={}, rec={}", orig.0, rec.0);
75        }
76    }
77
78    #[test]
79    fn spectral_heat_smooths() {
80        let n = 32;
81        let mut u: Vec<f64> = (0..n).map(|i| if i == n / 2 { 1.0 } else { 0.0 }).collect();
82        let peak_before = u[n / 2];
83        spectral_heat_step(&mut u, 1.0, 0.01, 1.0 / n as f64);
84        assert!(u[n / 2] < peak_before, "Heat should smooth the spike");
85    }
86}
proof_engine/solver/spectral.rs

proof_engine/solver/
spectral.rs
