diff --git a/russell_lab/data/figures/test_new_with_uu.svg b/russell_lab/data/figures/test_new_with_uu.svg index f0b821c3..d91ccc50 100644 --- a/russell_lab/data/figures/test_new_with_uu.svg +++ b/russell_lab/data/figures/test_new_with_uu.svg @@ -6,7 +6,7 @@ - 2024-06-25T17:48:19.778178 + 2024-06-26T11:48:51.565868 image/svg+xml @@ -42,16 +42,16 @@ z +" clip-path="url(#pa95b1753bb)" style="fill: none; stroke-dasharray: 2.96,1.28; stroke-dashoffset: 0; stroke: #808080; stroke-width: 0.8"/> - - + @@ -97,11 +97,11 @@ z +" clip-path="url(#pa95b1753bb)" style="fill: none; stroke-dasharray: 2.96,1.28; stroke-dashoffset: 0; stroke: #808080; stroke-width: 0.8"/> - + @@ -143,11 +143,11 @@ z +" clip-path="url(#pa95b1753bb)" style="fill: none; stroke-dasharray: 2.96,1.28; stroke-dashoffset: 0; stroke: #808080; stroke-width: 0.8"/> - + @@ -184,11 +184,11 @@ z +" clip-path="url(#pa95b1753bb)" style="fill: none; stroke-dasharray: 2.96,1.28; stroke-dashoffset: 0; stroke: #808080; stroke-width: 0.8"/> - + @@ -236,11 +236,11 @@ z +" clip-path="url(#pa95b1753bb)" style="fill: none; stroke-dasharray: 2.96,1.28; stroke-dashoffset: 0; stroke: #808080; stroke-width: 0.8"/> - + @@ -297,11 +297,11 @@ z +" clip-path="url(#pa95b1753bb)" style="fill: none; stroke-dasharray: 2.96,1.28; stroke-dashoffset: 0; stroke: #808080; stroke-width: 0.8"/> - + @@ -358,16 +358,16 @@ z +" clip-path="url(#pa95b1753bb)" style="fill: none; stroke-dasharray: 2.96,1.28; stroke-dashoffset: 0; stroke: #808080; stroke-width: 0.8"/> - - + @@ -381,11 +381,11 @@ L -3.5 0 +" clip-path="url(#pa95b1753bb)" style="fill: none; stroke-dasharray: 2.96,1.28; stroke-dashoffset: 0; stroke: #808080; stroke-width: 0.8"/> - + @@ -400,11 +400,11 @@ L 397.723125 233.816365 +" clip-path="url(#pa95b1753bb)" style="fill: none; stroke-dasharray: 2.96,1.28; stroke-dashoffset: 0; stroke: #808080; stroke-width: 0.8"/> - + @@ -419,11 +419,11 @@ L 397.723125 183.416365 +" clip-path="url(#pa95b1753bb)" style="fill: none; stroke-dasharray: 2.96,1.28; stroke-dashoffset: 0; stroke: #808080; stroke-width: 0.8"/> - + @@ -472,11 +472,11 @@ z +" clip-path="url(#pa95b1753bb)" style="fill: none; stroke-dasharray: 2.96,1.28; stroke-dashoffset: 0; stroke: #808080; stroke-width: 0.8"/> - + @@ -491,11 +491,11 @@ L 397.723125 82.616365 +" clip-path="url(#pa95b1753bb)" style="fill: none; stroke-dasharray: 2.96,1.28; stroke-dashoffset: 0; stroke: #808080; stroke-width: 0.8"/> - + @@ -634,51 +634,54 @@ L 344.155125 56.923845 L 363.634398 49.180897 L 381.490398 42.296365 L 381.490398 42.296365 -" clip-path="url(#pd3de1e6443)" style="fill: none; stroke: #1f77b4; stroke-width: 1.5; stroke-linecap: square"/> +" clip-path="url(#pa95b1753bb)" style="fill: none; stroke: #1f77b4; stroke-width: 1.5; stroke-linecap: square"/> +" clip-path="url(#pa95b1753bb)" style="fill: none; stroke-dasharray: 1.5,2.475; stroke-dashoffset: 0; stroke: #ff7f0e; stroke-width: 1.5"/> - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +" clip-path="url(#pa95b1753bb)" style="fill: none; stroke: #808080; stroke-width: 1.5; stroke-linecap: square"/> +" clip-path="url(#pa95b1753bb)" style="fill: none; stroke: #808080; stroke-width: 1.5; stroke-linecap: square"/> - - + - - + + - - + - + diff --git a/russell_lab/src/algo/interp_chebyshev.rs b/russell_lab/src/algo/interp_chebyshev.rs index 1b748fe5..f8205d29 100644 --- a/russell_lab/src/algo/interp_chebyshev.rs +++ b/russell_lab/src/algo/interp_chebyshev.rs @@ -1,4 +1,4 @@ -use crate::math::{chebyshev_lobatto_points_standard, chebyshev_tn, PI}; +use crate::math::{chebyshev_lobatto_points, chebyshev_tn, PI}; use crate::StrError; use crate::{NoArgs, Vector}; @@ -41,7 +41,6 @@ pub const TOL_RANGE: f64 = 1.0e-5; /// the derivative matrices, whereas this structure does not. /// 2. The [crate::InterpLagrange] renders the same results when using Chebyshev-Gauss-Lobatto points. /// 3. Only Chebyshev-Gauss-Lobatto points are considered here. -/// 4. The Chebyshev-Gauss-Lobatto coordinates are sorted from +1 to -1 (as in Reference # 1). /// /// # References /// @@ -63,11 +62,15 @@ pub struct InterpChebyshev { /// Holds the difference xb - xa dx: f64, - /// Holds the expansion coefficients (standard Chebyshev-Gauss-Lobatto) + /// Holds the expansion coefficients + /// + /// (associated with the reversed (from 1 to -1) Chebyshev-Gauss-Lobatto points) a: Vector, - /// Holds the function evaluation at the standard (from 1 to -1) Chebyshev-Gauss-Lobatto points - uu: Vector, + /// Holds the reversed function evaluations + /// + /// (associated with the reversed (from 1 to -1) Chebyshev-Gauss-Lobatto points) + uu_rev: Vector, /// Holds the constant y=c value for a zeroth-order function constant_fx: f64, @@ -77,13 +80,13 @@ pub struct InterpChebyshev { } impl InterpChebyshev { - /// Returns the standard (from 1 to -1) Chebyshev-Gauss-Lobatto points + /// Returns the Chebyshev-Gauss-Lobatto points (from -1 to 1) /// /// # Input /// /// * `nn` -- polynomial degree N pub fn points(nn: usize) -> Vector { - chebyshev_lobatto_points_standard(nn) + chebyshev_lobatto_points(nn) } /// Allocates a new instance with uninitialized values @@ -108,7 +111,7 @@ impl InterpChebyshev { xb, dx: xb - xa, a: Vector::new(np), - uu: Vector::new(np), + uu_rev: Vector::new(np), constant_fx: 0.0, ready: false, }) @@ -153,9 +156,9 @@ impl InterpChebyshev { /// } /// ``` pub fn set_uu_value(&mut self, i: usize, uui: f64) { - self.uu[i] = uui; + self.uu_rev[self.nn - i] = uui; if i == self.nn { - chebyshev_coefficients(self.a.as_mut_data(), self.uu.as_mut_data(), self.nn); + chebyshev_coefficients(self.a.as_mut_data(), self.uu_rev.as_mut_data(), self.nn); self.ready = true; } else { self.ready = false; @@ -207,15 +210,15 @@ impl InterpChebyshev { xb, dx: xb - xa, a: Vector::new(np), - uu: Vector::new(np), + uu_rev: Vector::new(np), constant_fx: 0.0, ready: true, }; if nn == 0 { interp.constant_fx = f((xa + xb) / 2.0, args)?; } else { - chebyshev_data_vector(interp.uu.as_mut_data(), nn, xa, xb, args, &mut f)?; - chebyshev_coefficients(interp.a.as_mut_data(), interp.uu.as_mut_data(), nn); + chebyshev_data_vector(interp.uu_rev.as_mut_data(), nn, xa, xb, args, &mut f)?; + chebyshev_coefficients(interp.a.as_mut_data(), interp.uu_rev.as_mut_data(), nn); } Ok(interp) } @@ -227,7 +230,7 @@ impl InterpChebyshev { /// * `xa` -- lower bound /// * `xb` -- upper bound (> xa + ϵ) /// * `uu` -- the data vector such that `Uᵢ = f(xᵢ)`; i.e., the function evaluated at the - /// **standard** (from 1 to -1) Chebyshev-Gauss-Lobatto coordinates. These coordinates + /// Chebyshev-Gauss-Lobatto coordinates (from -1 to 1). These coordinates /// are available via the [InterpChebyshev::points()] function. pub fn new_with_uu(xa: f64, xb: f64, uu: &[f64]) -> Result { if xb <= xa + TOL_RANGE { @@ -245,14 +248,17 @@ impl InterpChebyshev { xb, dx: xb - xa, a: Vector::new(np), - uu: Vector::from(&uu), + uu_rev: Vector::new(np), constant_fx: 0.0, ready: true, }; + for i in 0..np { + interp.uu_rev[nn - i] = uu[i]; + } if nn == 0 { interp.constant_fx = uu[0]; } else { - chebyshev_coefficients(interp.a.as_mut_data(), interp.uu.as_mut_data(), nn); + chebyshev_coefficients(interp.a.as_mut_data(), interp.uu_rev.as_mut_data(), nn); } Ok(interp) } @@ -320,11 +326,11 @@ impl InterpChebyshev { } let np_max = nn_max + 1; let mut work_a = vec![0.0; np_max]; - let mut work_uu = vec![0.0; np_max]; + let mut work_uu_rev = vec![0.0; np_max]; let mut an_prev = 0.0; for nn in 1..=nn_max { - chebyshev_data_vector(&mut work_uu, nn, xa, xb, args, &mut f)?; - chebyshev_coefficients(&mut work_a, &work_uu, nn); + chebyshev_data_vector(&mut work_uu_rev, nn, xa, xb, args, &mut f)?; + chebyshev_coefficients(&mut work_a, &work_uu_rev, nn); let an = work_a[nn]; if nn > 1 && f64::abs(an_prev) < tol && f64::abs(an) < tol { let nn_final = nn - 2; // -2 because the last two coefficients are zero @@ -343,7 +349,9 @@ impl InterpChebyshev { /// * `tol` -- tolerance to truncate the Chebyshev series (e.g., 1e-8) /// * `xa` -- lower bound /// * `xb` -- upper bound (> xa + ϵ) - /// * `uu` -- the data vector (len > 0) + /// * `uu` -- the data vector such that `Uᵢ = f(xᵢ)`; i.e., the function evaluated at the + /// Chebyshev-Gauss-Lobatto coordinates (from -1 to 1). These coordinates + /// are available via the [InterpChebyshev::points()] function. /// /// # Method /// @@ -356,7 +364,7 @@ impl InterpChebyshev { /// /// fn main() -> Result<(), StrError> { /// // data - /// let uu = [3.0, 0.5, -4.5, -7.0]; + /// let uu = [-7.0, -4.0, 0.5, 3.0]; /// let (xa, xb) = (0.0, 1.0); /// /// // interpolant @@ -471,7 +479,7 @@ impl InterpChebyshev { /// Computes the data vector (function evaluations at Chebyshev-Gauss-Lobatto points) fn chebyshev_data_vector( - work_uu: &mut [f64], + work_uu_rev: &mut [f64], nn: usize, xa: f64, xb: f64, @@ -485,31 +493,31 @@ where let np = nn + 1; assert!(nn > 0); assert!(xb > xa); - assert!(work_uu.len() >= np); + assert!(work_uu_rev.len() >= np); - // data vector U + // reverse data vector U (associated to Chebyshev-Gauss-Lobatto points from 1 to -1) let nf = nn as f64; - let uu = &mut work_uu[0..np]; + let uu_rev = &mut work_uu_rev[0..np]; for k in 0..np { let kf = k as f64; let x = (xb + xa + (xb - xa) * f64::cos(PI * kf / nf)) / 2.0; - uu[k] = (*f)(x, args)?; + uu_rev[k] = (*f)(x, args)?; } Ok(()) } /// Computes the Chebyshev-Gauss-Lobatto coefficients -fn chebyshev_coefficients(work_a: &mut [f64], work_uu: &[f64], nn: usize) { +fn chebyshev_coefficients(work_a: &mut [f64], work_uu_rev: &[f64], nn: usize) { // check let np = nn + 1; assert!(nn > 0); assert!(work_a.len() >= np); - assert!(work_uu.len() >= np); + assert!(work_uu_rev.len() >= np); // coefficients a let nf = nn as f64; let a = &mut work_a[0..np]; - let uu = &work_uu[0..np]; + let uu_rev = &work_uu_rev[0..np]; for j in 0..np { let jf = j as f64; let qj = if j == 0 || j == nn { 2.0 } else { 1.0 }; @@ -517,7 +525,7 @@ fn chebyshev_coefficients(work_a: &mut [f64], work_uu: &[f64], nn: usize) { for k in 0..np { let kf = k as f64; let qk = if k == 0 || k == nn { 2.0 } else { 1.0 }; - a[j] += uu[k] * 2.0 * f64::cos(PI * jf * kf / nf) / (qj * qk * nf); + a[j] += uu_rev[k] * 2.0 * f64::cos(PI * jf * kf / nf) / (qj * qk * nf); } } } @@ -553,7 +561,7 @@ mod tests { assert_eq!(interp.nn, nn); assert_eq!(interp.np, np); assert_eq!(interp.a.dim(), np); - assert_eq!(interp.uu.dim(), np); + assert_eq!(interp.uu_rev.dim(), np); assert_eq!(interp.constant_fx, 0.0); assert_eq!(interp.ready, false); } @@ -617,7 +625,7 @@ mod tests { let f = |x, _: &mut NoArgs| Ok(-1.0 + f64::sqrt(1.0 + 2.0 * x * 1200.0)); let (xa, xb) = (0.0, 1.0); - let nn = 10; + let nn = 6; let zz = InterpChebyshev::points(nn); let dx = xb - xa; let args = &mut 0; @@ -632,11 +640,11 @@ mod tests { for i in 0..np { let x = (xb + xa + dx * zz[i]) / 2.0; let fxi = interp.eval(x).unwrap(); - approx_eq(fxi, interp.uu[i], 1e-13); + approx_eq(fxi, interp.uu_rev[nn - i], 1e-13); } let err = interp.estimate_max_error(100, args, f).unwrap(); println!("err = {}", err); - assert!(err < 0.73); + approx_eq(err, 1.74, 1e-3); // plot f(x) /* diff --git a/russell_lab/src/algo/multi_root_solver_cheby.rs b/russell_lab/src/algo/multi_root_solver_cheby.rs index 229c7309..9fd07d78 100644 --- a/russell_lab/src/algo/multi_root_solver_cheby.rs +++ b/russell_lab/src/algo/multi_root_solver_cheby.rs @@ -616,7 +616,7 @@ mod tests { fn linear_function_no_roots_works() { // data let (xa, xb) = (0.0, 1.0); - let uu = Vector::from(&[3.0, 0.5]); + let uu = Vector::from(&[0.5, 3.0]); // interpolant let nn_max = 100; @@ -636,7 +636,7 @@ mod tests { // data let (xa, xb) = (0.0, 1.0); let dx = xb - xa; - let uu = Vector::from(&[3.0, 0.5, -4.5, -7.0]); + let uu = Vector::from(&[-7.0, -4.5, 0.5, 3.0]); let np = uu.dim(); // number of points let nn = np - 1; // degree let mut xx_dat = Vector::new(np); diff --git a/russell_lab/src/math/chebyshev.rs b/russell_lab/src/math/chebyshev.rs index 38e21805..7c151e0d 100644 --- a/russell_lab/src/math/chebyshev.rs +++ b/russell_lab/src/math/chebyshev.rs @@ -330,114 +330,11 @@ pub fn chebyshev_lobatto_points(nn: usize) -> Vector { xx } -/// Returns the standard (from 1 to -1) Chebyshev-Gauss-Lobatto coordinates -/// -/// The point coordinates are defined by: -/// -/// ```text -/// ⎛ j⋅π ⎞ -/// Xⱼ = cos ⎜ ————— ⎟ -/// ⎝ N ⎠ -/// -/// j = 0 ... N -/// ``` -/// -/// Nonetheless, here, the coordinates are calculates by using the sin(x) function: -/// -/// ```text -/// ⎛ π⋅(N-2j) ⎞ -/// Xⱼ = sin ⎜ —————————— ⎟ -/// ⎝ 2⋅N ⎠ -/// -/// j = 0 ... N -/// ``` -/// -/// See equation (2.4.14) on page 86 of the Reference #1 (with the "standard" ordering of nodes from +1 to -1). -/// See also equation (1) of the Reference #2 -/// -/// **Note:** This function enforces the symmetry of the sequence of points. -/// -/// # Input -/// -/// * `nn` -- the polynomial degree `N`; thus the number of points is `npoint = nn + 1`. -/// -/// # Input -/// -/// * The number of points will be equal to `npoint = yy.len()` and the -/// polynomial degree will be equal to `nn = npoint - 1` -/// -/// # Output -/// -/// * `yy` -- is the array holding the coordinates (from 1 to -1). Its length is equal -/// to `npoint = nn + 1` where `nn` is the polynomial degree. Requirement: `npoint ≥ 2`. -/// -/// # Notes -/// -/// See the note below from Reference # 1 (page 86): -/// -/// "Note that the Chebyshev quadrature points as just defined are ordered -/// from right to left. This violates our general convention that quadrature points -/// are ordered from left to right (see Sect. 2.2.3). Virtually all of the classical -/// literature on Chebyshev spectral methods uses this reversed order. Therefore, -/// in the special case of the Chebyshev quadrature points we shall adhere to the -/// ordering convention that is widely used in the literature (and implemented -/// in the available software). We realize that our resolution of this dilemma -/// imposes upon the reader the task of mentally reversing the ordering of the -/// Chebyshev nodes whenever they are used in general formulas for orthogonal -/// polynomials." (Canuto, Hussaini, Quarteroni, Zang) -/// -/// # References -/// -/// 1. Canuto C, Hussaini MY, Quarteroni A, Zang TA (2006) Spectral Methods: Fundamentals in -/// Single Domains. Springer. 563p -/// 2. Baltensperger R and Trummer MR (2003) Spectral differencing with a twist, -/// SIAM Journal Scientific Computation 24(5):1465-1487 -/// -/// # Examples -/// -/// ![Chebyshev points](https://raw.githubusercontent.com/cpmech/russell/main/russell_lab/data/figures/math_chebyshev_points.svg) -/// -/// ``` -/// use russell_lab::math; -/// -/// let xx = math::chebyshev_lobatto_points_standard(8); -/// println!("\nChebyshev-Gauss-Lobatto points =\n{:.3?}\n", xx.as_data()); -/// ``` -/// -/// The output looks like: -/// -/// ```text -/// Chebyshev-Gauss-Lobatto points = -/// [1.000, 0.924, 0.707, 0.383, 0.000, -0.383, -0.707, -0.924, -1.000] -/// ``` -pub fn chebyshev_lobatto_points_standard(nn: usize) -> Vector { - let mut xx = Vector::new(nn + 1); - xx[0] = 1.0; - xx[nn] = -1.0; - if nn < 3 { - return xx; - } - let nf = nn as f64; - let d = 2.0 * nf; - let l = if (nn & 1) == 0 { - // even number of segments - nn / 2 - } else { - // odd number of segments - (nn + 3) / 2 - 1 - }; - for i in 1..l { - xx[nn - i] = -f64::sin(PI * (nf - 2.0 * (i as f64)) / d); - xx[i] = -xx[nn - i]; - } - return xx; -} - //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// #[cfg(test)] mod tests { - use super::{chebyshev_gauss_points, chebyshev_lobatto_points, chebyshev_lobatto_points_standard}; + use super::{chebyshev_gauss_points, chebyshev_lobatto_points}; use super::{chebyshev_tn, chebyshev_tn_deriv1, chebyshev_tn_deriv2}; use crate::math::{SQRT_2, SQRT_3}; use crate::{approx_eq, vec_approx_eq, Vector}; @@ -740,106 +637,4 @@ mod tests { vec_approx_eq(&xx, &xx_ref, 1e-15); check_segment_symmetry(&xx); } - - #[test] - fn chebyshev_lobatto_points_standard_works() { - let xx = chebyshev_lobatto_points_standard(0); - assert_eq!(xx.as_data(), &[-1.0]); - - let xx = chebyshev_lobatto_points_standard(1); - let xx_ref = vec![1.0, -1.0]; - vec_approx_eq(&xx, &xx_ref, 1e-15); - check_segment_symmetry(&xx); - - let xx = chebyshev_lobatto_points_standard(2); - let xx_ref = vec![1.0, 0.0, -1.0]; - vec_approx_eq(&xx, &xx_ref, 1e-15); - check_segment_symmetry(&xx); - - let xx = chebyshev_lobatto_points_standard(3); - let xx_ref = vec![1.0, 0.5, -0.5, -1.0]; - vec_approx_eq(&xx, &xx_ref, 1e-15); - check_segment_symmetry(&xx); - - let xx = chebyshev_lobatto_points_standard(4); - let xx_ref = vec![ - 1.000000000000000e+00, - 7.071067811865476e-01, - 0.0, - -7.071067811865475e-01, - -1.000000000000000e+00, - ]; - vec_approx_eq(&xx, &xx_ref, 1e-15); - check_segment_symmetry(&xx); - - let xx = chebyshev_lobatto_points_standard(5); - let xx_ref = vec![ - 1.000000000000000e+00, - 8.090169943749475e-01, - 3.090169943749475e-01, - -3.090169943749473e-01, - -8.090169943749473e-01, - -1.000000000000000e+00, - ]; - vec_approx_eq(&xx, &xx_ref, 1e-15); - check_segment_symmetry(&xx); - - let xx = chebyshev_lobatto_points_standard(6); - let xx_ref = vec![ - 1.000000000000000e+00, - 8.660254037844387e-01, - 5.000000000000001e-01, - 0.0, - -4.999999999999998e-01, - -8.660254037844385e-01, - -1.000000000000000e+00, - ]; - vec_approx_eq(&xx, &xx_ref, 1e-15); - check_segment_symmetry(&xx); - - let xx = chebyshev_lobatto_points_standard(7); - let xx_ref = vec![ - 1.000000000000000e+00, - 9.009688679024191e-01, - 6.234898018587336e-01, - 2.225209339563144e-01, - -2.225209339563143e-01, - -6.234898018587335e-01, - -9.009688679024190e-01, - -1.000000000000000e+00, - ]; - vec_approx_eq(&xx, &xx_ref, 1e-15); - check_segment_symmetry(&xx); - - let xx = chebyshev_lobatto_points_standard(8); - let xx_ref = vec![ - 1.000000000000000e+00, - 9.238795325112867e-01, - 7.071067811865476e-01, - 3.826834323650898e-01, - 0.0, - -3.826834323650897e-01, - -7.071067811865475e-01, - -9.238795325112867e-01, - -1.000000000000000e+00, - ]; - vec_approx_eq(&xx, &xx_ref, 1e-15); - check_segment_symmetry(&xx); - - let xx = chebyshev_lobatto_points_standard(9); - let xx_ref = vec![ - 1.000000000000000e+00, - 9.396926207859084e-01, - 7.660444431189780e-01, - 5.000000000000001e-01, - 1.736481776669304e-01, - -1.736481776669303e-01, - -4.999999999999998e-01, - -7.660444431189779e-01, - -9.396926207859083e-01, - -1.000000000000000e+00, - ]; - vec_approx_eq(&xx, &xx_ref, 1e-15); - check_segment_symmetry(&xx); - } }