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
pub fn chebyquad(x: &[f64], m: usize) -> f64 {
let n = x.len();
if m < n {
panic!("number of auxiliary function must be at least n");
}
if n == 0 {
panic!("input vector x must not be empty");
}
let n_f64 = n as f64;
let mut f_results = Vec::with_capacity(m);
// Loop i from 1 to m (mathematical index) to calculate each f_i(x)
for i in 1..=m {
// --- Calculate Summation Part ---
// (1/n) * Σ_{j=1 to n} T_i(x_j)
let mut sum_part = 0.0;
for xj in x {
sum_part += shifted_chebyshev_t(i, *xj);
}
sum_part /= n_f64;
// --- Calculate Integral Part ---
let integral_part = integral_ti(i);
// --- Combine and store the result for f_i ---
f_results.push(sum_part - integral_part);
}
f_results.iter().map(|&f_i| f_i.powi(2)).sum()
}
pub fn init(n: usize) -> Vec<f64> {
let mut vec = vec![0.; n];
for (i, item) in vec.iter_mut().enumerate().take(n) {
*item = (i as f64 + 1.) / (n as f64 + 1.0)
}
vec
}
/// Calculates the value of the i-th Chebyshev polynomial of the first kind,
/// shifted to the interval [0, 1].
/// Standard Chebyshev polynomials T_i(z) are defined on [-1, 1].
/// To shift them to x in [0, 1], we use the transformation z = 2x - 1.
///
/// They follow the recurrence relation:
/// T_0(z) = 1
/// T_1(z) = z
/// T_i(z) = 2z * T_{i-1}(z) - T_{i-2}(z)
fn shifted_chebyshev_t(i: usize, x: f64) -> f64 {
// Map x from [0, 1] to z in [-1, 1]
let z = 2.0 * x - 1.0;
if i == 0 {
return 1.0;
}
if i == 1 {
return z;
}
// Use the recurrence relation for i >= 2
let mut t_prev = 1.0; // T_0(z)
let mut t_curr = z; // T_1(z)
for _ in 2..=i {
let t_next = 2.0 * z * t_curr - t_prev;
t_prev = t_curr;
t_curr = t_next;
}
t_curr
}
/// Calculates the definite integral of the i-th shifted Chebyshev polynomial
/// from 0 to 1, based on the pre-calculated formula in the image.
fn integral_ti(i: usize) -> f64 {
// For i=1 (odd), we correctly return 0. The even case i^2-1 would be a division by zero.
if !i.is_multiple_of(2) {
// for i odd
0.0
} else {
// for i even
-1.0 / ((i as f64).powi(2) - 1.0)
}
}
#[cfg(test)]
mod tests {
use super::*;
#[test]
fn test_chebyquad() {
let n = 4;
let x = init(n);
let val = chebyquad(&x, 4);
assert!(val.is_finite());
}
}