use oxilean_kernel::{BinderInfo, Declaration, Environment, Expr, Level, Name};
use super::types::{
BoundedOp, Distribution, FredholmOperatorData, InterpolationData, InterpolationMethod,
L2Sequence, LuDecomposition, QrDecomposition, RnVector, SobolevSpaceData, WeakConvergenceData,
};
pub fn app(f: Expr, a: Expr) -> Expr {
Expr::App(Box::new(f), Box::new(a))
}
pub fn app2(f: Expr, a: Expr, b: Expr) -> Expr {
app(app(f, a), b)
}
pub fn cst(s: &str) -> Expr {
Expr::Const(Name::str(s), vec![])
}
pub fn prop() -> Expr {
Expr::Sort(Level::zero())
}
pub fn type0() -> Expr {
Expr::Sort(Level::succ(Level::zero()))
}
pub fn pi(bi: BinderInfo, name: &str, dom: Expr, body: Expr) -> Expr {
Expr::Pi(bi, Name::str(name), Box::new(dom), Box::new(body))
}
pub fn arrow(a: Expr, b: Expr) -> Expr {
pi(BinderInfo::Default, "_", a, b)
}
pub fn bvar(n: u32) -> Expr {
Expr::BVar(n)
}
pub fn nat_ty() -> Expr {
cst("Nat")
}
pub fn real_ty() -> Expr {
cst("Real")
}
pub fn impl_pi(name: &str, dom: Expr, body: Expr) -> Expr {
pi(BinderInfo::Implicit, name, dom, body)
}
pub fn _use_helpers() {
let _ = (
app(cst("a"), cst("b")),
app2(cst("a"), cst("b"), cst("c")),
prop(),
type0(),
arrow(prop(), prop()),
bvar(0),
nat_ty(),
real_ty(),
impl_pi("x", type0(), bvar(0)),
);
}
pub fn normed_space_ty() -> Expr {
arrow(type0(), prop())
}
pub fn banach_space_ty() -> Expr {
arrow(type0(), prop())
}
pub fn hilbert_space_ty() -> Expr {
arrow(type0(), prop())
}
pub fn bounded_linear_op_ty() -> Expr {
pi(
BinderInfo::Default,
"X",
type0(),
pi(BinderInfo::Default, "Y", type0(), type0()),
)
}
pub fn dual_space_ty() -> Expr {
arrow(type0(), type0())
}
pub fn spectrum_ty() -> Expr {
pi(
BinderInfo::Default,
"X",
type0(),
arrow(
app2(cst("BoundedLinearOp"), bvar(0), bvar(0)),
app(cst("Set"), real_ty()),
),
)
}
pub fn fredholm_operator_ty() -> Expr {
pi(
BinderInfo::Default,
"X",
type0(),
pi(BinderInfo::Default, "Y", type0(), prop()),
)
}
pub fn fredholm_index_ty() -> Expr {
pi(
BinderInfo::Default,
"X",
type0(),
pi(
BinderInfo::Default,
"Y",
type0(),
arrow(app2(cst("BoundedLinearOp"), bvar(1), bvar(0)), cst("Int")),
),
)
}
pub fn compact_operator_ty() -> Expr {
pi(
BinderInfo::Default,
"X",
type0(),
pi(
BinderInfo::Default,
"Y",
type0(),
arrow(app2(cst("BoundedLinearOp"), bvar(1), bvar(0)), prop()),
),
)
}
pub fn hahn_banach_ty() -> Expr {
pi(
BinderInfo::Default,
"X",
type0(),
arrow(app(cst("NormedSpace"), bvar(0)), prop()),
)
}
pub fn open_mapping_ty() -> Expr {
pi(
BinderInfo::Default,
"X",
type0(),
pi(
BinderInfo::Default,
"Y",
type0(),
arrow(
app2(cst("BoundedLinearOp"), bvar(1), bvar(0)),
arrow(
app(cst("BanachSpace"), bvar(2)),
arrow(app(cst("BanachSpace"), bvar(1)), prop()),
),
),
),
)
}
pub fn closed_graph_ty() -> Expr {
pi(
BinderInfo::Default,
"X",
type0(),
pi(
BinderInfo::Default,
"Y",
type0(),
arrow(
arrow(bvar(1), bvar(0)),
arrow(app(cst("BanachSpace"), bvar(2)), prop()),
),
),
)
}
pub fn uniform_boundedness_ty() -> Expr {
pi(
BinderInfo::Default,
"X",
type0(),
pi(
BinderInfo::Default,
"Y",
type0(),
arrow(
app(cst("BanachSpace"), bvar(1)),
arrow(app(cst("BanachSpace"), bvar(0)), prop()),
),
),
)
}
pub fn riesz_representation_ty() -> Expr {
pi(
BinderInfo::Default,
"H",
type0(),
arrow(app(cst("HilbertSpace"), bvar(0)), prop()),
)
}
pub fn spectral_theorem_ty() -> Expr {
pi(
BinderInfo::Default,
"H",
type0(),
arrow(
app(cst("HilbertSpace"), bvar(0)),
arrow(app2(cst("CompactSelfAdjoint"), bvar(1), bvar(0)), prop()),
),
)
}
pub fn fredholm_alternative_ty() -> Expr {
pi(
BinderInfo::Default,
"X",
type0(),
arrow(
app(cst("BanachSpace"), bvar(0)),
arrow(app2(cst("BoundedLinearOp"), bvar(1), bvar(1)), prop()),
),
)
}
pub fn banach_fixed_point_ty() -> Expr {
pi(
BinderInfo::Default,
"X",
type0(),
arrow(
app(cst("BanachSpace"), bvar(0)),
arrow(arrow(bvar(1), bvar(1)), prop()),
),
)
}
pub fn build_functional_analysis_env(env: &mut Environment) {
let axioms: &[(&str, Expr)] = &[
("NormedSpace", normed_space_ty()),
("BanachSpace", banach_space_ty()),
("HilbertSpace", hilbert_space_ty()),
("BoundedLinearOp", bounded_linear_op_ty()),
("DualSpace", dual_space_ty()),
("Spectrum", spectrum_ty()),
("FredholmOperator", fredholm_operator_ty()),
("FredholmIndex", fredholm_index_ty()),
("CompactOperator", compact_operator_ty()),
("Set", arrow(type0(), type0())),
(
"CompactSelfAdjoint",
pi(BinderInfo::Default, "H", type0(), arrow(type0(), type0())),
),
("Int", type0()),
("hahn_banach", hahn_banach_ty()),
("open_mapping", open_mapping_ty()),
("closed_graph", closed_graph_ty()),
("uniform_boundedness", uniform_boundedness_ty()),
("riesz_representation", riesz_representation_ty()),
("spectral_theorem", spectral_theorem_ty()),
("fredholm_alternative", fredholm_alternative_ty()),
("banach_fixed_point", banach_fixed_point_ty()),
];
for (name, ty) in axioms {
env.add(Declaration::Axiom {
name: Name::str(*name),
univ_params: vec![],
ty: ty.clone(),
})
.ok();
}
}
pub fn gram_schmidt(vectors: &[RnVector]) -> Vec<RnVector> {
let mut orthonormal = Vec::new();
for v in vectors {
let mut u = v.clone();
for q in &orthonormal {
let proj = u.project_onto(q);
u = u.sub(&proj);
}
let n = u.norm();
if n > 1e-12 {
orthonormal.push(u.scale(1.0 / n));
}
}
orthonormal
}
pub fn is_orthonormal(vectors: &[RnVector], tol: f64) -> bool {
for (i, vi) in vectors.iter().enumerate() {
if (vi.norm() - 1.0).abs() > tol {
return false;
}
for vj in vectors.iter().skip(i + 1) {
if vi.inner(vj).abs() > tol {
return false;
}
}
}
true
}
pub fn spans_full_space(vectors: &[RnVector], dim: usize) -> bool {
gram_schmidt(vectors).len() == dim
}
pub fn conjugate_gradient(a: &BoundedOp, b: &RnVector, max_iter: usize, tol: f64) -> RnVector {
let n = b.dim();
let mut x = RnVector::zero(n);
let mut r = b.sub(&a.apply(&x));
let mut p = r.clone();
let mut rs_old = r.inner(&r);
for _ in 0..max_iter {
if rs_old.sqrt() < tol {
break;
}
let ap = a.apply(&p);
let alpha = rs_old / p.inner(&ap);
x = x.add(&p.scale(alpha));
r = r.sub(&ap.scale(alpha));
let rs_new = r.inner(&r);
if rs_new.sqrt() < tol {
break;
}
p = r.add(&p.scale(rs_new / rs_old));
rs_old = rs_new;
}
x
}
pub fn l2_norm_approx(f: &dyn Fn(f64) -> f64, a: f64, b: f64, n: usize) -> f64 {
if n == 0 {
return 0.0;
}
let h = (b - a) / n as f64;
let mut sum = 0.5 * f(a).powi(2) + 0.5 * f(b).powi(2);
for i in 1..n {
sum += f(a + i as f64 * h).powi(2);
}
(sum * h).sqrt()
}
pub fn l2_inner_approx(
f: &dyn Fn(f64) -> f64,
g: &dyn Fn(f64) -> f64,
a: f64,
b: f64,
n: usize,
) -> f64 {
if n == 0 {
return 0.0;
}
let h = (b - a) / n as f64;
let mut sum = 0.5 * f(a) * g(a) + 0.5 * f(b) * g(b);
for i in 1..n {
let x = a + i as f64 * h;
sum += f(x) * g(x);
}
sum * h
}
pub fn sup_norm_approx(f: &dyn Fn(f64) -> f64, a: f64, b: f64, n: usize) -> f64 {
if n == 0 {
return f(a).abs();
}
let h = (b - a) / n as f64;
let mut max_val = 0.0_f64;
for i in 0..=n {
max_val = max_val.max(f(a + i as f64 * h).abs());
}
max_val
}
pub fn banach_fixed_point_iter(
t: &dyn Fn(f64) -> f64,
x0: f64,
max_iter: usize,
tol: f64,
) -> (f64, usize) {
let mut x = x0;
for i in 0..max_iter {
let x_next = t(x);
if (x_next - x).abs() < tol {
return (x_next, i + 1);
}
x = x_next;
}
(x, max_iter)
}
pub fn fourier_cosine_coefficients(
f: &dyn Fn(f64) -> f64,
period: f64,
num_coeffs: usize,
qp: usize,
) -> Vec<f64> {
let h = period / qp as f64;
let mut coeffs = Vec::with_capacity(num_coeffs);
for k in 0..num_coeffs {
let mut sum =
0.5 * f(0.0) + 0.5 * f(period) * (2.0 * std::f64::consts::PI * k as f64).cos();
for i in 1..qp {
let x = i as f64 * h;
sum += f(x) * (2.0 * std::f64::consts::PI * k as f64 * x / period).cos();
}
let factor = if k == 0 { 1.0 / period } else { 2.0 / period };
coeffs.push(sum * h * factor);
}
coeffs
}
pub fn eval_fourier_cosine(coeffs: &[f64], period: f64, x: f64) -> f64 {
coeffs
.iter()
.enumerate()
.map(|(k, &c)| c * (2.0 * std::f64::consts::PI * k as f64 * x / period).cos())
.sum()
}
pub fn qr_decomposition(a: &BoundedOp) -> Option<QrDecomposition> {
let (m, n) = (a.range_dim, a.domain_dim);
if m == 0 || n == 0 {
return None;
}
let mut cols: Vec<RnVector> = (0..n)
.map(|j| RnVector::new((0..m).map(|i| a.matrix[i][j]).collect()))
.collect();
let mut q_cols: Vec<RnVector> = Vec::new();
let mut r_matrix = vec![vec![0.0; n]; n.min(m)];
for j in 0..n.min(m) {
for k in 0..q_cols.len() {
let rk = cols[j].inner(&q_cols[k]);
r_matrix[k][j] = rk;
cols[j] = cols[j].sub(&q_cols[k].scale(rk));
}
let nrm = cols[j].norm();
if nrm < 1e-14 {
r_matrix[j][j] = 0.0;
q_cols.push(RnVector::zero(m));
} else {
r_matrix[j][j] = nrm;
q_cols.push(cols[j].scale(1.0 / nrm));
}
}
let q_n = q_cols.len();
let q_matrix: Vec<Vec<f64>> = (0..m)
.map(|i| (0..q_n).map(|j| q_cols[j].components[i]).collect())
.collect();
let q = BoundedOp {
matrix: q_matrix,
domain_dim: q_n,
range_dim: m,
};
let r = BoundedOp {
matrix: r_matrix,
domain_dim: n,
range_dim: q_n,
};
Some(QrDecomposition { q, r })
}
pub fn lu_decomposition(a: &BoundedOp) -> Option<LuDecomposition> {
let n = a.domain_dim;
if n != a.range_dim || n == 0 {
return None;
}
let mut u = a.matrix.clone();
let mut l = vec![vec![0.0; n]; n];
let mut perm: Vec<usize> = (0..n).collect();
for k in 0..n {
let mut max_val = u[k][k].abs();
let mut max_row = k;
for row in (k + 1)..n {
if u[row][k].abs() > max_val {
max_val = u[row][k].abs();
max_row = row;
}
}
if max_val < 1e-15 {
return None;
}
if max_row != k {
u.swap(k, max_row);
l.swap(k, max_row);
perm.swap(k, max_row);
}
l[k][k] = 1.0;
for i in (k + 1)..n {
let factor = u[i][k] / u[k][k];
l[i][k] = factor;
for j in k..n {
let val = u[k][j];
u[i][j] -= factor * val;
}
}
}
Some(LuDecomposition {
l: BoundedOp {
matrix: l,
domain_dim: n,
range_dim: n,
},
u: BoundedOp {
matrix: u,
domain_dim: n,
range_dim: n,
},
permutation: perm,
})
}
pub fn lu_solve(lu: &LuDecomposition, b: &RnVector) -> RnVector {
let n = lu.l.domain_dim;
let mut pb = vec![0.0; n];
for (i, &pi) in lu.permutation.iter().enumerate() {
pb[i] = b.components[pi];
}
let mut y = vec![0.0; n];
for i in 0..n {
let mut sum = pb[i];
for j in 0..i {
sum -= lu.l.matrix[i][j] * y[j];
}
y[i] = sum;
}
let mut x = vec![0.0; n];
for i in (0..n).rev() {
let mut sum = y[i];
for j in (i + 1)..n {
sum -= lu.u.matrix[i][j] * x[j];
}
if lu.u.matrix[i][i].abs() < 1e-15 {
x[i] = 0.0;
} else {
x[i] = sum / lu.u.matrix[i][i];
}
}
RnVector::new(x)
}
pub fn singular_values(a: &BoundedOp, max_iter: usize) -> Vec<f64> {
let ata = match a.transpose().compose(a) {
Some(m) => m,
None => return vec![],
};
let n = ata.domain_dim;
if n == 0 {
return vec![];
}
let mut current = ata.clone();
for _ in 0..max_iter {
let qr = match qr_decomposition(¤t) {
Some(qr) => qr,
None => break,
};
current = match qr.r.compose(&qr.q) {
Some(m) => m,
None => break,
};
}
let mut eigenvalues: Vec<f64> = (0..n).map(|i| current.matrix[i][i].max(0.0)).collect();
eigenvalues.sort_by(|a, b| b.partial_cmp(a).unwrap_or(std::cmp::Ordering::Equal));
eigenvalues.iter().map(|&ev| ev.sqrt()).collect()
}
pub fn condition_number(a: &BoundedOp, max_iter: usize) -> f64 {
let svs = singular_values(a, max_iter);
if svs.is_empty() {
return f64::INFINITY;
}
let max_sv = svs[0];
let min_sv = *svs
.last()
.expect("svs is non-empty: checked by early return");
if min_sv < 1e-15 {
f64::INFINITY
} else {
max_sv / min_sv
}
}
pub fn spectral_radius(a: &BoundedOp, iterations: usize) -> f64 {
if a.domain_dim != a.range_dim || a.domain_dim == 0 {
return 0.0;
}
let n = a.domain_dim;
let mut v = RnVector::new(vec![1.0; n]);
let nrm = v.norm();
if nrm > 0.0 {
v = v.scale(1.0 / nrm);
}
let mut radius = 0.0;
for _ in 0..iterations {
let w = a.apply(&v);
let nrm = w.norm();
if nrm < 1e-15 {
return 0.0;
}
radius = nrm;
v = w.scale(1.0 / nrm);
}
radius
}
pub fn h1_norm_approx(f: &dyn Fn(f64) -> f64, a: f64, b: f64, n: usize) -> f64 {
let l2_sq = l2_norm_approx(f, a, b, n).powi(2);
if n < 2 {
return l2_sq.sqrt();
}
let h = (b - a) / n as f64;
let mut deriv_sq_sum = 0.0;
for i in 0..n {
let x = a + i as f64 * h;
let df = (f(x + h) - f(x)) / h;
deriv_sq_sum += df * df;
}
let deriv_l2_sq = deriv_sq_sum * h;
(l2_sq + deriv_l2_sq).sqrt()
}
pub fn h2_norm_approx(f: &dyn Fn(f64) -> f64, a: f64, b: f64, n: usize) -> f64 {
let h1_sq = h1_norm_approx(f, a, b, n).powi(2);
if n < 3 {
return h1_sq.sqrt();
}
let h = (b - a) / n as f64;
let mut d2_sq_sum = 0.0;
for i in 1..(n - 1) {
let x = a + i as f64 * h;
let d2f = (f(x + h) - 2.0 * f(x) + f(x - h)) / (h * h);
d2_sq_sum += d2f * d2f;
}
let d2_l2_sq = d2_sq_sum * h;
(h1_sq + d2_l2_sq).sqrt()
}
pub fn simpson_integrate(f: &dyn Fn(f64) -> f64, a: f64, b: f64, n: usize) -> f64 {
let n = if n % 2 == 0 { n } else { n + 1 };
if n == 0 {
return 0.0;
}
let h = (b - a) / n as f64;
let mut sum = f(a) + f(b);
for i in 1..n {
let x = a + i as f64 * h;
sum += if i % 2 == 0 { 2.0 * f(x) } else { 4.0 * f(x) };
}
sum * h / 3.0
}
pub fn gauss_legendre_3(f: &dyn Fn(f64) -> f64, a: f64, b: f64) -> f64 {
let mid = (a + b) / 2.0;
let half = (b - a) / 2.0;
let nodes = [-0.7745966692414834, 0.0, 0.7745966692414834];
let weights = [0.5555555555555556, 0.8888888888888888, 0.5555555555555556];
let mut sum = 0.0;
for i in 0..3 {
sum += weights[i] * f(mid + half * nodes[i]);
}
sum * half
}
pub fn gauss_legendre_composite(f: &dyn Fn(f64) -> f64, a: f64, b: f64, panels: usize) -> f64 {
if panels == 0 {
return 0.0;
}
let h = (b - a) / panels as f64;
let mut total = 0.0;
for i in 0..panels {
let lo = a + i as f64 * h;
let hi = lo + h;
total += gauss_legendre_3(f, lo, hi);
}
total
}
pub fn matrix_exponential(a: &BoundedOp, terms: usize) -> BoundedOp {
let n = a.domain_dim;
if n != a.range_dim {
return BoundedOp::zero_op(n, n);
}
let mut result = BoundedOp::identity(n);
let mut term = BoundedOp::identity(n);
for k in 1..terms {
term = match term.compose(a) {
Some(m) => m.scalar_mul(1.0 / k as f64),
None => break,
};
result = match result.op_add(&term) {
Some(m) => m,
None => break,
};
}
result
}
pub fn commutator(a: &BoundedOp, b: &BoundedOp) -> Option<BoundedOp> {
let ab = a.compose(b)?;
let ba = b.compose(a)?;
let neg_ba = ba.scalar_mul(-1.0);
ab.op_add(&neg_ba)
}
pub fn anti_commutator(a: &BoundedOp, b: &BoundedOp) -> Option<BoundedOp> {
let ab = a.compose(b)?;
let ba = b.compose(a)?;
ab.op_add(&ba)
}
pub fn fredholm_index_numerical(k: &BoundedOp) -> (usize, usize, i64) {
let n = k.domain_dim;
if n != k.range_dim {
return (0, 0, 0);
}
let i_minus_k = match BoundedOp::identity(n).op_add(&k.scalar_mul(-1.0)) {
Some(m) => m,
None => return (0, 0, 0),
};
let kernel_dim = i_minus_k.nullity();
let cokernel_dim = i_minus_k.transpose().nullity();
let index = kernel_dim as i64 - cokernel_dim as i64;
(kernel_dim, cokernel_dim, index)
}
pub fn jacobi_iteration(a: &BoundedOp, b: &RnVector, max_iter: usize, tol: f64) -> RnVector {
let n = a.domain_dim;
if n != a.range_dim || n != b.dim() {
return RnVector::zero(n);
}
let mut x = RnVector::zero(n);
for _ in 0..max_iter {
let mut x_new = vec![0.0; n];
for i in 0..n {
if a.matrix[i][i].abs() < 1e-15 {
x_new[i] = x.components[i];
continue;
}
let mut sum = b.components[i];
for j in 0..n {
if j != i {
sum -= a.matrix[i][j] * x.components[j];
}
}
x_new[i] = sum / a.matrix[i][i];
}
let new_x = RnVector::new(x_new);
if new_x.sub(&x).norm() < tol {
return new_x;
}
x = new_x;
}
x
}
pub fn gauss_seidel_iteration(a: &BoundedOp, b: &RnVector, max_iter: usize, tol: f64) -> RnVector {
let n = a.domain_dim;
if n != a.range_dim || n != b.dim() {
return RnVector::zero(n);
}
let mut x = vec![0.0; n];
for _ in 0..max_iter {
let old_x = x.clone();
for i in 0..n {
if a.matrix[i][i].abs() < 1e-15 {
continue;
}
let mut sum = b.components[i];
for j in 0..n {
if j != i {
sum -= a.matrix[i][j] * x[j];
}
}
x[i] = sum / a.matrix[i][i];
}
let diff: f64 = x
.iter()
.zip(old_x.iter())
.map(|(a, b)| (a - b).powi(2))
.sum::<f64>()
.sqrt();
if diff < tol {
break;
}
}
RnVector::new(x)
}
pub fn fourier_sine_coefficients(
f: &dyn Fn(f64) -> f64,
period: f64,
num_coeffs: usize,
qp: usize,
) -> Vec<f64> {
let h = period / qp as f64;
let mut coeffs = Vec::with_capacity(num_coeffs);
for k in 0..num_coeffs {
let mut sum = 0.0;
for i in 0..=qp {
let x = i as f64 * h;
let w = if i == 0 || i == qp { 0.5 } else { 1.0 };
sum += w * f(x) * (std::f64::consts::PI * (k + 1) as f64 * x / period).sin();
}
coeffs.push(sum * h * 2.0 / period);
}
coeffs
}
pub fn eval_fourier_sine(coeffs: &[f64], period: f64, x: f64) -> f64 {
coeffs
.iter()
.enumerate()
.map(|(k, &c)| c * (std::f64::consts::PI * (k + 1) as f64 * x / period).sin())
.sum()
}
pub fn chebyshev_t(n: usize, x: f64) -> f64 {
if n == 0 {
return 1.0;
}
if n == 1 {
return x;
}
let (mut t0, mut t1) = (1.0, x);
for _ in 2..=n {
let t2 = 2.0 * x * t1 - t0;
t0 = t1;
t1 = t2;
}
t1
}
pub fn chebyshev_nodes(n: usize, a: f64, b: f64) -> Vec<f64> {
(0..n)
.map(|k| {
let theta = std::f64::consts::PI * (2 * k + 1) as f64 / (2 * n) as f64;
(a + b) / 2.0 + (b - a) / 2.0 * theta.cos()
})
.collect()
}
pub fn chebyshev_coefficients(f: &dyn Fn(f64) -> f64, n: usize, a: f64, b: f64) -> Vec<f64> {
let nodes = chebyshev_nodes(n, a, b);
let values: Vec<f64> = nodes.iter().map(|&x| f(x)).collect();
let mut coeffs = vec![0.0; n];
for j in 0..n {
let mut sum = 0.0;
for (k, &val) in values.iter().enumerate() {
let xk = 2.0 * (nodes[k] - a) / (b - a) - 1.0;
sum += val * chebyshev_t(j, xk);
}
coeffs[j] = sum
* if j == 0 {
1.0 / n as f64
} else {
2.0 / n as f64
};
}
coeffs
}
pub fn eval_chebyshev(coeffs: &[f64], a: f64, b: f64, x: f64) -> f64 {
let t = 2.0 * (x - a) / (b - a) - 1.0;
coeffs
.iter()
.enumerate()
.map(|(k, &c)| c * chebyshev_t(k, t))
.sum()
}
#[cfg(test)]
mod tests {
use super::*;
#[test]
fn test_rn_vector_norm() {
let v = RnVector::new(vec![3.0, 4.0]);
assert!((v.norm() - 5.0).abs() < 1e-10);
}
#[test]
fn test_rn_vector_inner() {
let e1 = RnVector::new(vec![1.0, 0.0]);
let e2 = RnVector::new(vec![0.0, 1.0]);
assert!(e1.inner(&e2).abs() < 1e-10);
assert!((e1.inner(&e1) - 1.0).abs() < 1e-10);
}
#[test]
fn test_rn_vector_normalized() {
let v = RnVector::new(vec![3.0, 4.0]);
let u = v.normalized().expect("nonzero");
assert!((u.norm() - 1.0).abs() < 1e-10);
assert!(RnVector::zero(3).normalized().is_none());
}
#[test]
fn test_rn_vector_basis() {
let e0 = RnVector::basis(3, 0);
assert!((e0.components[0] - 1.0).abs() < 1e-10);
assert!(e0.components[1].abs() < 1e-10);
}
#[test]
fn test_rn_vector_sub() {
let c = RnVector::new(vec![5.0, 3.0]).sub(&RnVector::new(vec![2.0, 1.0]));
assert!((c.components[0] - 3.0).abs() < 1e-10);
assert!((c.components[1] - 2.0).abs() < 1e-10);
}
#[test]
fn test_rn_vector_cross() {
let k = RnVector::new(vec![1.0, 0.0, 0.0])
.cross(&RnVector::new(vec![0.0, 1.0, 0.0]))
.expect("cross should succeed");
assert!(k.components[0].abs() < 1e-10);
assert!(k.components[1].abs() < 1e-10);
assert!((k.components[2] - 1.0).abs() < 1e-10);
}
#[test]
fn test_rn_vector_project() {
let proj = RnVector::new(vec![3.0, 4.0]).project_onto(&RnVector::new(vec![1.0, 0.0]));
assert!((proj.components[0] - 3.0).abs() < 1e-10);
assert!(proj.components[1].abs() < 1e-10);
}
#[test]
fn test_gram_schmidt() {
let basis = gram_schmidt(&[
RnVector::new(vec![1.0, 1.0, 0.0]),
RnVector::new(vec![1.0, 0.0, 1.0]),
RnVector::new(vec![0.0, 1.0, 1.0]),
]);
assert_eq!(basis.len(), 3);
assert!(is_orthonormal(&basis, 1e-10));
}
#[test]
fn test_gram_schmidt_dependent() {
let basis = gram_schmidt(&[
RnVector::new(vec![1.0, 0.0]),
RnVector::new(vec![2.0, 0.0]),
RnVector::new(vec![0.0, 1.0]),
]);
assert_eq!(basis.len(), 2);
assert!(is_orthonormal(&basis, 1e-10));
}
#[test]
fn test_spans_full_space() {
assert!(spans_full_space(
&[RnVector::new(vec![1.0, 0.0]), RnVector::new(vec![0.0, 1.0])],
2
));
assert!(!spans_full_space(
&[
RnVector::new(vec![1.0, 0.0, 0.0]),
RnVector::new(vec![0.0, 1.0, 0.0])
],
3
));
}
#[test]
fn test_bounded_op_apply() {
let id = BoundedOp::identity(3);
let v = RnVector::new(vec![1.0, 2.0, 3.0]);
let w = id.apply(&v);
for (a, b) in w.components.iter().zip(v.components.iter()) {
assert!((a - b).abs() < 1e-10);
}
}
#[test]
fn test_bounded_op_transpose() {
let a = BoundedOp::new(vec![vec![1.0, 2.0, 3.0], vec![4.0, 5.0, 6.0]]);
let at = a.transpose();
assert_eq!(at.range_dim, 3);
assert_eq!(at.domain_dim, 2);
let att = at.transpose();
for i in 0..2 {
for j in 0..3 {
assert!((att.matrix[i][j] - a.matrix[i][j]).abs() < 1e-10);
}
}
}
#[test]
fn test_bounded_op_symmetric() {
assert!(BoundedOp::new(vec![vec![2.0, 1.0], vec![1.0, 3.0]]).is_symmetric());
assert!(!BoundedOp::new(vec![vec![1.0, 2.0], vec![3.0, 4.0]]).is_symmetric());
}
#[test]
fn test_bounded_op_diagonal() {
let w = BoundedOp::diagonal(&[2.0, 3.0, 5.0]).apply(&RnVector::new(vec![1.0, 1.0, 1.0]));
assert!((w.components[0] - 2.0).abs() < 1e-10);
assert!((w.components[1] - 3.0).abs() < 1e-10);
assert!((w.components[2] - 5.0).abs() < 1e-10);
}
#[test]
fn test_bounded_op_add() {
let c = BoundedOp::identity(2)
.op_add(&BoundedOp::identity(2))
.expect("op_add should succeed");
assert!((c.matrix[0][0] - 2.0).abs() < 1e-10);
}
#[test]
fn test_bounded_op_determinant() {
assert!(
(BoundedOp::identity(3)
.determinant()
.expect("determinant should succeed")
- 1.0)
.abs()
< 1e-10
);
assert!(
(BoundedOp::new(vec![vec![1.0, 2.0], vec![3.0, 4.0]])
.determinant()
.expect("determinant should succeed")
- (-2.0))
.abs()
< 1e-10
);
}
#[test]
fn test_bounded_op_solve() {
let x = BoundedOp::new(vec![vec![2.0, 1.0], vec![1.0, 3.0]])
.solve(&RnVector::new(vec![5.0, 7.0]))
.expect("solve should succeed");
assert!((x.components[0] - 1.6).abs() < 1e-10);
assert!((x.components[1] - 1.8).abs() < 1e-10);
}
#[test]
fn test_bounded_op_rank_nullity() {
let a = BoundedOp::new(vec![
vec![1.0, 2.0, 3.0],
vec![4.0, 5.0, 6.0],
vec![7.0, 8.0, 9.0],
]);
assert_eq!(a.rank(), 2);
assert_eq!(a.nullity(), 1);
}
#[test]
fn test_bounded_op_positive_definite() {
assert!(BoundedOp::new(vec![vec![2.0, 1.0], vec![1.0, 2.0]]).is_positive_definite());
assert!(!BoundedOp::new(vec![vec![1.0, 2.0], vec![2.0, 1.0]]).is_positive_definite());
}
#[test]
fn test_power_iteration() {
let (ev, _) = BoundedOp::diagonal(&[5.0, 1.0])
.power_iteration(100)
.expect("power_iteration should succeed");
assert!((ev - 5.0).abs() < 1e-6, "got {ev}");
}
#[test]
fn test_eigenvalues_2x2() {
let (e1, e2) = BoundedOp::new(vec![vec![3.0, 1.0], vec![1.0, 3.0]])
.eigenvalues_2x2()
.expect("eigenvalues_2x2 should succeed");
assert!((e1 - 4.0).abs() < 1e-10);
assert!((e2 - 2.0).abs() < 1e-10);
}
#[test]
fn test_conjugate_gradient() {
let a = BoundedOp::new(vec![vec![4.0, 1.0], vec![1.0, 3.0]]);
let x = conjugate_gradient(&a, &RnVector::new(vec![1.0, 2.0]), 100, 1e-12);
let ax = a.apply(&x);
assert!((ax.components[0] - 1.0).abs() < 1e-8);
assert!((ax.components[1] - 2.0).abs() < 1e-8);
}
#[test]
fn test_conjugate_gradient_identity() {
let b = RnVector::new(vec![1.0, 2.0, 3.0]);
let x = conjugate_gradient(&BoundedOp::identity(3), &b, 100, 1e-12);
for i in 0..3 {
assert!((x.components[i] - b.components[i]).abs() < 1e-8);
}
}
#[test]
fn test_l2_sequence_norm() {
assert!((L2Sequence::new(vec![1.0, 0.0, 0.0]).l2_norm() - 1.0).abs() < 1e-10);
let b = L2Sequence::new(vec![3.0, 4.0]);
assert!((b.l2_norm() - 5.0).abs() < 1e-10);
assert!(b.is_in_l2());
}
#[test]
fn test_l2_sequence_shift() {
let s = L2Sequence::new(vec![1.0, 2.0, 3.0]);
assert_eq!(s.shift_left().terms, vec![2.0, 3.0]);
assert_eq!(s.shift_right().terms, vec![0.0, 1.0, 2.0, 3.0]);
}
#[test]
fn test_l2_sequence_convolve() {
let c = L2Sequence::new(vec![1.0, 2.0]).convolve(&L2Sequence::new(vec![3.0, 4.0]));
assert!((c.terms[0] - 3.0).abs() < 1e-10);
assert!((c.terms[1] - 10.0).abs() < 1e-10);
assert!((c.terms[2] - 8.0).abs() < 1e-10);
}
#[test]
fn test_l2_sequence_parseval() {
assert!(L2Sequence::new(vec![1.0, 2.0, 3.0, 4.0]).parseval_residual() < 1e-10);
}
#[test]
fn test_l2_norm_approx() {
assert!((l2_norm_approx(&|_| 1.0, 0.0, 1.0, 1000) - 1.0).abs() < 1e-3);
}
#[test]
fn test_l2_inner_approx() {
let ip = l2_inner_approx(
&|x: f64| x.sin(),
&|x: f64| x.cos(),
0.0,
2.0 * std::f64::consts::PI,
10000,
);
assert!(ip.abs() < 1e-3, "got {ip}");
}
#[test]
fn test_sup_norm_approx() {
let n = sup_norm_approx(&|x: f64| x.sin(), 0.0, std::f64::consts::PI, 10000);
assert!((n - 1.0).abs() < 1e-3, "got {n}");
}
#[test]
fn test_banach_fixed_point_iter() {
let (fp, iters) = banach_fixed_point_iter(&|x: f64| x.cos(), 0.5, 1000, 1e-10);
assert!((fp - fp.cos()).abs() < 1e-8);
assert!(iters < 1000);
}
#[test]
fn test_fourier_cosine_coefficients() {
let c = fourier_cosine_coefficients(&|_| 1.0, 1.0, 5, 1000);
assert!((c[0] - 1.0).abs() < 1e-3);
for k in 1..5 {
assert!(c[k].abs() < 1e-3);
}
}
#[test]
fn test_eval_fourier_cosine() {
assert!((eval_fourier_cosine(&[1.0, 0.0, 0.0], 1.0, 0.5) - 1.0).abs() < 1e-10);
}
#[test]
fn test_bounded_op_operator_norm_power() {
let n = BoundedOp::diagonal(&[3.0, 1.0, 2.0]).operator_norm_power_iter(200);
assert!((n - 3.0).abs() < 1e-3, "got {n}");
}
#[test]
fn test_bounded_op_frobenius_norm() {
let frob = BoundedOp::identity(2).frobenius_norm();
assert!((frob - 2.0_f64.sqrt()).abs() < 1e-10);
assert!(BoundedOp::identity(2).operator_norm() > 0.0);
}
#[test]
fn test_bounded_op_compose() {
let a = BoundedOp::new(vec![vec![1.0, 2.0], vec![3.0, 4.0]]);
let c = a
.compose(&BoundedOp::identity(2))
.expect("BoundedOp::new should succeed");
for i in 0..2 {
for j in 0..2 {
assert!((c.matrix[i][j] - a.matrix[i][j]).abs() < 1e-10);
}
}
}
#[test]
fn test_bounded_op_scalar_mul() {
let b = BoundedOp::identity(2).scalar_mul(3.0);
assert!((b.matrix[0][0] - 3.0).abs() < 1e-10);
}
#[test]
fn test_rn_vector_lp_norms() {
let v = RnVector::new(vec![3.0, 4.0]);
assert!((v.norm_p(1.0) - 7.0).abs() < 1e-10);
assert!((v.norm_p(f64::INFINITY) - 4.0).abs() < 1e-10);
}
#[test]
fn test_build_functional_analysis_env() {
let mut env = Environment::new();
build_functional_analysis_env(&mut env);
assert!(env.get(&Name::str("NormedSpace")).is_some());
assert!(env.get(&Name::str("hahn_banach")).is_some());
assert!(env.get(&Name::str("fredholm_alternative")).is_some());
assert!(env.get(&Name::str("banach_fixed_point")).is_some());
}
#[test]
fn test_bounded_op_injective_surjective() {
assert!(BoundedOp::identity(3).is_injective());
assert!(BoundedOp::identity(3).is_surjective());
let s = BoundedOp::new(vec![vec![1.0, 0.0], vec![0.0, 0.0]]);
assert!(!s.is_injective());
assert!(!s.is_surjective());
}
#[test]
fn test_rn_vector_angle() {
let angle = RnVector::new(vec![1.0, 0.0]).angle_with(&RnVector::new(vec![0.0, 1.0]));
assert!((angle - std::f64::consts::FRAC_PI_2).abs() < 1e-10);
}
#[test]
fn test_qr_identity() {
let qr = qr_decomposition(&BoundedOp::identity(3)).expect("operation should succeed");
for i in 0..3 {
for j in 0..3 {
let expected = if i == j { 1.0 } else { 0.0 };
assert!((qr.q.matrix[i][j] - expected).abs() < 1e-10);
assert!((qr.r.matrix[i][j] - expected).abs() < 1e-10);
}
}
}
#[test]
fn test_qr_reconstruction() {
let a = BoundedOp::new(vec![vec![1.0, 2.0], vec![3.0, 4.0], vec![5.0, 6.0]]);
let qr = qr_decomposition(&a).expect("BoundedOp::new should succeed");
let prod = qr.q.compose(&qr.r).expect("compose should succeed");
for i in 0..a.range_dim {
for j in 0..a.domain_dim {
assert!(
(prod.matrix[i][j] - a.matrix[i][j]).abs() < 1e-10,
"mismatch at ({i},{j}): {} vs {}",
prod.matrix[i][j],
a.matrix[i][j]
);
}
}
}
#[test]
fn test_qr_orthogonality() {
let a = BoundedOp::new(vec![
vec![12.0, -51.0, 4.0],
vec![6.0, 167.0, -68.0],
vec![-4.0, 24.0, -41.0],
]);
let qr = qr_decomposition(&a).expect("operation should succeed");
let qtq =
qr.q.transpose()
.compose(&qr.q)
.expect("compose should succeed");
for i in 0..3 {
for j in 0..3 {
let expected = if i == j { 1.0 } else { 0.0 };
assert!((qtq.matrix[i][j] - expected).abs() < 1e-10);
}
}
}
#[test]
fn test_lu_solve() {
let a = BoundedOp::new(vec![
vec![2.0, 1.0, 1.0],
vec![4.0, 3.0, 3.0],
vec![8.0, 7.0, 9.0],
]);
let b = RnVector::new(vec![1.0, 1.0, 1.0]);
let lu = lu_decomposition(&a).expect("RnVector::new should succeed");
let x = lu_solve(&lu, &b);
let ax = a.apply(&x);
for i in 0..3 {
assert!(
(ax.components[i] - b.components[i]).abs() < 1e-8,
"lu_solve error at {i}"
);
}
}
#[test]
fn test_lu_identity() {
let lu = lu_decomposition(&BoundedOp::identity(3)).expect("operation should succeed");
for i in 0..3 {
assert!((lu.l.matrix[i][i] - 1.0).abs() < 1e-10);
assert!((lu.u.matrix[i][i] - 1.0).abs() < 1e-10);
}
}
#[test]
fn test_singular_values_diagonal() {
let svs = singular_values(&BoundedOp::diagonal(&[3.0, 1.0, 5.0]), 50);
assert!((svs[0] - 5.0).abs() < 0.1, "sv[0] = {}", svs[0]);
assert!((svs[1] - 3.0).abs() < 0.1, "sv[1] = {}", svs[1]);
assert!((svs[2] - 1.0).abs() < 0.1, "sv[2] = {}", svs[2]);
}
#[test]
fn test_condition_number_identity() {
let cn = condition_number(&BoundedOp::identity(3), 50);
assert!((cn - 1.0).abs() < 0.1, "cond(I) = {cn}");
}
#[test]
fn test_condition_number_ill_conditioned() {
let cn = condition_number(&BoundedOp::diagonal(&[100.0, 0.01]), 50);
assert!(cn > 1000.0, "expected ill-conditioned, got {cn}");
}
#[test]
fn test_spectral_radius() {
let rho = spectral_radius(&BoundedOp::diagonal(&[2.0, -3.0, 1.0]), 200);
assert!((rho - 3.0).abs() < 0.1, "got {rho}");
}
#[test]
fn test_h1_norm_constant() {
let h1 = h1_norm_approx(&|_| 1.0, 0.0, 1.0, 1000);
assert!((h1 - 1.0).abs() < 1e-2, "got {h1}");
}
#[test]
fn test_h1_norm_linear() {
let h1 = h1_norm_approx(&|x| x, 0.0, 1.0, 1000);
let expected = (1.0 / 3.0 + 1.0_f64).sqrt();
assert!(
(h1 - expected).abs() < 0.05,
"got {h1}, expected {expected}"
);
}
#[test]
fn test_h2_norm() {
let h2 = h2_norm_approx(&|x: f64| x * x, 0.0, 1.0, 1000);
assert!(h2 > 0.0);
}
#[test]
fn test_simpson_integrate() {
let result = simpson_integrate(&|x| x * x, 0.0, 1.0, 100);
assert!((result - 1.0 / 3.0).abs() < 1e-6, "got {result}");
}
#[test]
fn test_gauss_legendre_quadrature() {
let result = gauss_legendre_composite(&|x| x * x, 0.0, 1.0, 10);
assert!((result - 1.0 / 3.0).abs() < 1e-10, "got {result}");
}
#[test]
fn test_matrix_exponential_zero() {
let exp_zero = matrix_exponential(&BoundedOp::zero_op(3, 3), 20);
for i in 0..3 {
for j in 0..3 {
let expected = if i == j { 1.0 } else { 0.0 };
assert!((exp_zero.matrix[i][j] - expected).abs() < 1e-10);
}
}
}
#[test]
fn test_matrix_exponential_identity() {
let exp_i = matrix_exponential(&BoundedOp::identity(2), 20);
let e = std::f64::consts::E;
assert!(
(exp_i.matrix[0][0] - e).abs() < 1e-6,
"got {}",
exp_i.matrix[0][0]
);
assert!((exp_i.matrix[1][1] - e).abs() < 1e-6);
assert!(exp_i.matrix[0][1].abs() < 1e-6);
}
#[test]
fn test_commutator_self() {
let a = BoundedOp::new(vec![vec![1.0, 2.0], vec![3.0, 4.0]]);
let c = commutator(&a, &a).expect("BoundedOp::new should succeed");
for i in 0..2 {
for j in 0..2 {
assert!(c.matrix[i][j].abs() < 1e-10);
}
}
}
#[test]
fn test_anti_commutator() {
let a = BoundedOp::identity(2);
let ac = anti_commutator(&a, &a).expect("BoundedOp::identity should succeed");
assert!((ac.matrix[0][0] - 2.0).abs() < 1e-10);
assert!((ac.matrix[1][1] - 2.0).abs() < 1e-10);
}
#[test]
fn test_fredholm_index_zero() {
let k = BoundedOp::zero_op(3, 3);
let (kd, cd, idx) = fredholm_index_numerical(&k);
assert_eq!(kd, 0);
assert_eq!(cd, 0);
assert_eq!(idx, 0);
}
#[test]
fn test_fredholm_index_projection() {
let k = BoundedOp::new(vec![vec![1.0, 0.0], vec![0.0, 0.0]]);
let (kd, _cd, _idx) = fredholm_index_numerical(&k);
assert_eq!(kd, 1);
}
#[test]
fn test_jacobi_iteration() {
let a = BoundedOp::new(vec![vec![4.0, 1.0], vec![1.0, 3.0]]);
let b = RnVector::new(vec![1.0, 2.0]);
let x = jacobi_iteration(&a, &b, 200, 1e-10);
let ax = a.apply(&x);
assert!((ax.components[0] - b.components[0]).abs() < 1e-6);
assert!((ax.components[1] - b.components[1]).abs() < 1e-6);
}
#[test]
fn test_gauss_seidel_iteration() {
let a = BoundedOp::new(vec![vec![4.0, 1.0], vec![1.0, 3.0]]);
let b = RnVector::new(vec![1.0, 2.0]);
let x = gauss_seidel_iteration(&a, &b, 200, 1e-10);
let ax = a.apply(&x);
assert!((ax.components[0] - b.components[0]).abs() < 1e-6);
assert!((ax.components[1] - b.components[1]).abs() < 1e-6);
}
#[test]
fn test_fourier_sine_coefficients() {
let coeffs =
fourier_sine_coefficients(&|x: f64| (std::f64::consts::PI * x).sin(), 1.0, 5, 1000);
assert!((coeffs[0] - 1.0).abs() < 1e-2, "got {}", coeffs[0]);
for k in 1..5 {
assert!(coeffs[k].abs() < 1e-2, "coeff[{k}] = {}", coeffs[k]);
}
}
#[test]
fn test_eval_fourier_sine() {
assert!((eval_fourier_sine(&[1.0, 0.0, 0.0], 1.0, 0.5) - 1.0).abs() < 1e-10);
}
#[test]
fn test_chebyshev_t_values() {
assert!((chebyshev_t(0, 0.5) - 1.0).abs() < 1e-10);
assert!((chebyshev_t(1, 0.5) - 0.5).abs() < 1e-10);
assert!((chebyshev_t(2, 0.5) - (-0.5)).abs() < 1e-10);
}
#[test]
fn test_chebyshev_nodes() {
let nodes = chebyshev_nodes(4, -1.0, 1.0);
assert_eq!(nodes.len(), 4);
for &x in &nodes {
assert!(x >= -1.0 && x <= 1.0);
}
}
#[test]
fn test_chebyshev_approximation() {
let coeffs = chebyshev_coefficients(&|x: f64| x.sin(), 10, -1.0, 1.0);
let approx = eval_chebyshev(&coeffs, -1.0, 1.0, 0.5);
assert!(
(approx - 0.5_f64.sin()).abs() < 0.1,
"got {approx}, expected {}",
0.5_f64.sin()
);
}
}
#[cfg(test)]
mod tests_functional_analysis_ext {
use super::*;
#[test]
fn test_fredholm_operator() {
let fo = FredholmOperatorData::new("D_A", 2, 1);
assert_eq!(fo.index(), 1);
assert!(!fo.is_isomorphism());
assert!(fo.index_stable_under_compact());
assert!(fo.atkinson_description().contains("index 1"));
}
#[test]
fn test_fredholm_isomorphism() {
let fo = FredholmOperatorData::new("I", 0, 0);
assert_eq!(fo.index(), 0);
assert!(fo.is_isomorphism());
}
#[test]
fn test_sobolev_space() {
let ss = SobolevSpaceData::new(1, 2.0, 3, "R^3");
let p_star = ss
.critical_sobolev_exponent()
.expect("critical_sobolev_exponent should succeed");
assert!((p_star - 6.0).abs() < 1e-10, "p* should be 6, got {p_star}");
assert!(ss.rellich_kondrachov_compact());
assert!(ss.trace_theorem().contains("Trace"));
}
#[test]
fn test_sobolev_supercritical() {
let ss = SobolevSpaceData::new(2, 2.0, 2, "R^2");
assert!(ss.critical_sobolev_exponent().is_none());
}
#[test]
fn test_interpolation() {
let id = InterpolationData::new("L^1", "L^inf", 0.5, InterpolationMethod::Complex);
let exp = id.lp_exponent(1.0, 1e10);
assert!((exp - 0.5).abs() < 1e-9, "1/p should be ~0.5");
let bound = id.riesz_thorin_bound(1.0, 1.0);
assert!((bound - 1.0).abs() < 1e-10);
}
}
#[cfg(test)]
mod tests_functional_analysis_ext2 {
use super::*;
#[test]
fn test_weak_convergence() {
let seq: Vec<Vec<f64>> = (1..=10).map(|n| vec![1.0 / n as f64, 0.0]).collect();
let wcd = WeakConvergenceData::new(seq);
assert!(wcd.is_bounded);
assert!(wcd.banach_alaoglu_applies());
assert!(wcd.check_weak_convergence(&[0.0, 0.0], 0.2));
}
#[test]
fn test_distribution_dirac() {
let delta = Distribution::dirac_delta(0.0);
assert!(!delta.is_regular);
assert_eq!(delta.order, Some(0));
assert!(delta.is_tempered());
let deriv = delta.differentiate();
assert_eq!(deriv.order, Some(1));
}
#[test]
fn test_distribution_regular() {
let reg = Distribution::regular("f", "R").differentiate();
assert!(!reg.is_regular);
assert_eq!(reg.order, Some(1));
assert!(reg
.fourier_transform_description()
.contains("distributional"));
}
}