From c0e70ef8953ee0f987a6b33b38e0baf52963ab06 Mon Sep 17 00:00:00 2001 From: Asfmq <2696428814@qq.com> Date: Wed, 25 Mar 2026 01:54:00 +0800 Subject: [PATCH] =?UTF-8?q?feat:=20=E6=B7=BB=E5=8A=A0=20lagran=20=E5=92=8C?= =?UTF-8?q?=20yint=20=E6=A8=A1=E5=9D=97?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - lagran: Lagrange 三点插值函数 - yint: 二次插值函数（数组版本） 这两个纯数学函数用于数值插值。 Co-Authored-By: Claude Opus 4.6 --- src/math/lagran.rs | 81 ++++++++++++++++++++++++++++++++++++++++++++++ src/math/mod.rs | 2 ++ src/math/yint.rs | 80 +++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 163 insertions(+) create mode 100644 src/math/lagran.rs create mode 100644 src/math/yint.rs diff --git a/src/math/lagran.rs b/src/math/lagran.rs new file mode 100644 index 0000000..4a05bcb --- /dev/null +++ b/src/math/lagran.rs @@ -0,0 +1,81 @@ +//! Lagrange 三点插值。 +//! +//! 重构自 TLUSTY `lagran.f` + +/// Lagrange 三点插值。 +/// +/// 给定三个点 (x0,y0), (x1,y1), (x2,y2)，计算 x 处的插值 y。 +/// +/// # 参数 +/// * `x0`, `y0` - 第一个点 +/// * `x1`, `y1` - 第二个点 +/// * `x2`, `y2` - 第三个点 +/// * `x` - 需要插值的位置 +/// +/// # 返回值 +/// 返回插值结果 y +/// +/// # 注意 +/// 调用者需确保 x0, x1, x2 互不相等 +pub fn lagran(x0: f64, x1: f64, x2: f64, y0: f64, y1: f64, y2: f64, x: f64) -> f64 { + // Lagrange 基函数 + let xl0 = (x - x1) * (x - x2) / ((x0 - x1) * (x0 - x2)); + let xl1 = (x - x0) * (x - x2) / ((x1 - x0) * (x1 - x2)); + let xl2 = (x - x0) * (x - x1) / ((x2 - x0) * (x2 - x1)); + + // 插值结果 + y0 * xl0 + y1 * xl1 + y2 * xl2 +} + +#[cfg(test)] +mod tests { + use super::*; + use approx::assert_relative_eq; + + #[test] + fn test_lagran_at_points() { + // 测试在插值节点上的值 + let y = lagran(0.0, 1.0, 2.0, 1.0, 2.0, 5.0, 0.0); + assert_relative_eq!(y, 1.0, epsilon = 1e-15); + + let y = lagran(0.0, 1.0, 2.0, 1.0, 2.0, 5.0, 1.0); + assert_relative_eq!(y, 2.0, epsilon = 1e-15); + + let y = lagran(0.0, 1.0, 2.0, 1.0, 2.0, 5.0, 2.0); + assert_relative_eq!(y, 5.0, epsilon = 1e-15); + } + + #[test] + fn test_lagran_interpolation() { + // 测试线性插值 + let y = lagran(0.0, 1.0, 2.0, 0.0, 1.0, 2.0, 0.5); + assert_relative_eq!(y, 0.5, epsilon = 1e-15); + + // 测试二次插值 (y = x^2) + let y = lagran(0.0, 1.0, 2.0, 0.0, 1.0, 4.0, 0.5); + assert_relative_eq!(y, 0.25, epsilon = 1e-15); + + let y = lagran(0.0, 1.0, 2.0, 0.0, 1.0, 4.0, 1.5); + assert_relative_eq!(y, 2.25, epsilon = 1e-15); + } + + #[test] + fn test_lagran_extrapolation() { + // 测试外推 + let y = lagran(0.0, 1.0, 2.0, 0.0, 1.0, 4.0, -1.0); + assert_relative_eq!(y, 1.0, epsilon = 1e-15); // (-1)^2 = 1 + + let y = lagran(0.0, 1.0, 2.0, 0.0, 1.0, 4.0, 3.0); + assert_relative_eq!(y, 9.0, epsilon = 1e-15); // 3^2 = 9 + } + + #[test] + fn test_lagran_negative_x() { + // 测试负 x 值 + let y = lagran(-2.0, 0.0, 2.0, 4.0, 0.0, 4.0, 0.0); + assert_relative_eq!(y, 0.0, epsilon = 1e-15); + + let y = lagran(-2.0, 0.0, 2.0, 4.0, 0.0, 4.0, -1.0); + assert_relative_eq!(y, 1.0, epsilon = 1e-15); + } +} diff --git a/src/math/mod.rs b/src/math/mod.rs index 826f5a0..16656ee 100644 --- a/src/math/mod.rs +++ b/src/math/mod.rs @@ -112,6 +112,7 @@ mod intlem; mod intxen; mod irc; mod interpolate; +mod lagran; mod laguer; mod lemini; mod levsol; @@ -267,6 +268,7 @@ mod voigte; mod wn; mod wnstor; mod xk2dop; +mod yint; mod ylintp; mod zmrho; diff --git a/src/math/yint.rs b/src/math/yint.rs new file mode 100644 index 0000000..fdcc094 --- /dev/null +++ b/src/math/yint.rs @@ -0,0 +1,80 @@ +//! 二次插值函数。 +//! +//! 重构自 TLUSTY `yint.f` + +/// 二次插值。 +/// +/// 给定三个点 (xl[0],yl[0]), (xl[1],yl[1]), (xl[2],yl[2])， +/// 计算 xl0 处的插值结果。 +/// +/// # 参数 +/// * `xl` - x 坐标数组（长度为 3） +/// * `yl` - y 坐标数组（长度为 3） +/// * `xl0` - 需要插值的位置 +/// +/// # 返回值 +/// 返回插值结果 +/// +/// # 注意 +/// 调用者需确保 xl[0], xl[1], xl[2] 互不相等 +pub fn yint(xl: &[f64], yl: &[f64], xl0: f64) -> f64 { + // Fortran 索引转换：XL(1) -> xl[0], XL(2) -> xl[1], XL(3) -> xl[2] + let a0 = (xl[1] - xl[0]) * (xl[2] - xl[1]) * (xl[2] - xl[0]); + let a1 = (xl0 - xl[1]) * (xl0 - xl[2]) * (xl[2] - xl[1]); + let a2 = (xl0 - xl[0]) * (xl[2] - xl0) * (xl[2] - xl[0]); + let a3 = (xl0 - xl[0]) * (xl0 - xl[1]) * (xl[1] - xl[0]); + + (yl[0] * a1 + yl[1] * a2 + yl[2] * a3) / a0 +} + +#[cfg(test)] +mod tests { + use super::*; + use approx::assert_relative_eq; + + #[test] + fn test_yint_at_points() { + let xl = [0.0, 1.0, 2.0]; + let yl = [1.0, 2.0, 5.0]; + + // 测试在插值节点上的值 + assert_relative_eq!(yint(&xl, &yl, 0.0), 1.0, epsilon = 1e-15); + assert_relative_eq!(yint(&xl, &yl, 1.0), 2.0, epsilon = 1e-15); + assert_relative_eq!(yint(&xl, &yl, 2.0), 5.0, epsilon = 1e-15); + } + + #[test] + fn test_yint_quadratic() { + // y = x^2 + let xl = [0.0, 1.0, 2.0]; + let yl = [0.0, 1.0, 4.0]; + + assert_relative_eq!(yint(&xl, &yl, 0.5), 0.25, epsilon = 1e-15); + assert_relative_eq!(yint(&xl, &yl, 1.5), 2.25, epsilon = 1e-15); + } + + #[test] + fn test_yint_extrapolation() { + let xl = [0.0, 1.0, 2.0]; + let yl = [0.0, 1.0, 4.0]; + + // 外推: y = x^2 + assert_relative_eq!(yint(&xl, &yl, -1.0), 1.0, epsilon = 1e-15); + assert_relative_eq!(yint(&xl, &yl, 3.0), 9.0, epsilon = 1e-15); + } + + #[test] + fn test_yint_vs_lagran() { + // yint 和 lagran 应该给出相同的结果 + use super::super::lagran::lagran; + + let xl = [0.0, 1.0, 2.0]; + let yl = [1.0, 3.0, 2.0]; + + for x in [0.0, 0.5, 1.0, 1.5, 2.0, 3.0].iter() { + let y1 = yint(&xl, &yl, *x); + let y2 = lagran(xl[0], xl[1], xl[2], yl[0], yl[1], yl[2], *x); + assert_relative_eq!(y1, y2, epsilon = 1e-15); + } + } +}