roots/numerical/
mod.rs

1// Copyright (c) 2015, Mikhail Vorotilov
2// All rights reserved.
3//
4// Redistribution and use in source and binary forms, with or without
5// modification, are permitted provided that the following conditions are met:
6//
7// * Redistributions of source code must retain the above copyright notice, this
8//   list of conditions and the following disclaimer.
9//
10// * Redistributions in binary form must reproduce the above copyright notice,
11//   this list of conditions and the following disclaimer in the documentation
12//   and/or other materials provided with the distribution.
13//
14// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
15// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
16// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
17// DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
18// FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
19// DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
20// SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
21// CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
22// OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
23// OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
24
25use super::FloatType;
26use std::error::Error;
27use std::fmt;
28
29/// Pair of the independent variable x and the function value y=F(x)
30#[derive(Debug, PartialEq)]
31pub struct Sample<F>
32where
33    F: FloatType,
34{
35    /// Value of the independent variable (X-axis)
36    x: F,
37    /// Value of the dependent variable (Y-axis)
38    y: F,
39}
40
41impl<F> Sample<F>
42where
43    F: FloatType,
44{
45    fn is_bracketed_with(&self, other: &Self) -> bool {
46        self.y * other.y <= F::zero()
47    }
48}
49
50/// Interval between two samples, including these samples
51#[derive(Debug, PartialEq)]
52pub struct Interval<F>
53where
54    F: FloatType,
55{
56    /// First sample
57    begin: Sample<F>,
58    /// Last sample
59    end: Sample<F>,
60}
61
62impl<F> Interval<F>
63where
64    F: FloatType,
65{
66    fn is_bracketed(&self) -> bool {
67        self.begin.is_bracketed_with(&self.end)
68    }
69    fn is_converged(&self, convergency: &mut dyn Convergency<F>) -> bool {
70        convergency.is_converged(self.begin.x, self.end.x)
71    }
72    /// Check if the given X is inside the interval
73    fn contains_x(&self, x: &F) -> bool {
74        *x <= self.end.x && *x >= self.begin.x
75    }
76    /// Returns a point somewhere in middle of the interval for narrowing this interval down.
77    /// Rules are as follows:
78    /// * If the interval is bracketed, use the secant to find the middle point.
79    /// ** The middle point may not be too close to either range of the interval.
80    /// * If the interval is not bracketed (why would one use an unbracketed interval?), bisect it.
81    fn middle(&self) -> F {
82        let _2 = F::from(2i16);
83        let _26 = F::from(26i16);
84        let _27 = F::from(27i16);
85
86        if self.is_bracketed() && self.begin.y != self.end.y {
87            let mut shift = -self.begin.y * (self.end.x - self.begin.x) / (self.end.y - self.begin.y);
88            if shift < (self.end.x - self.begin.x) / _27 {
89                shift = (self.end.x - self.begin.x) / _27;
90            }
91            if shift > (self.end.x - self.begin.x) * _26 / _27 {
92                shift = (self.end.x - self.begin.x) * _26 / _27;
93            }
94            self.begin.x + shift
95        } else {
96            (self.begin.x + self.end.x) / _2
97        }
98    }
99}
100
101/// Possible errors
102#[derive(Debug, PartialEq)]
103pub enum SearchError {
104    /// The algorithm could not converge within the given number of iterations
105    NoConvergency,
106    /// Initial values do not bracket zero
107    NoBracketing,
108    /// The algorithm cannot continue from the point where the derivative is zero
109    ZeroDerivative,
110}
111
112impl fmt::Display for SearchError {
113    fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
114        match self {
115            SearchError::NoConvergency => write!(f, "Convergency Error"),
116            SearchError::NoBracketing => write!(f, "Bracketing Error"),
117            SearchError::ZeroDerivative => write!(f, "Zero Derivative Error"),
118        }
119    }
120}
121impl Error for SearchError {
122    fn description(&self) -> &str {
123        match self {
124            SearchError::NoConvergency => "The algorithm could not converge within the given number of iterations",
125            SearchError::NoBracketing => "Initial values do not bracket zero",
126            SearchError::ZeroDerivative => "The algorithm cannot continue from the point where the derivative is zero",
127        }
128    }
129}
130
131/// The way to check if the algorithm has finished by either finding a root
132/// or reaching the iteration limit.
133pub trait Convergency<F: FloatType> {
134    /// Return true if the given Y value is close enough to the zero
135    fn is_root_found(&mut self, y: F) -> bool;
136    /// Return true if given x values are close enough to each other
137    fn is_converged(&mut self, x1: F, x2: F) -> bool;
138    /// Return true if no more iterations desired
139    fn is_iteration_limit_reached(&mut self, iter: usize) -> bool;
140}
141
142impl<F: FloatType> Convergency<F> for F {
143    /// Return true if the given Y value is close enough to the zero
144    fn is_root_found(&mut self, y: F) -> bool {
145        y.abs() < self.abs()
146    }
147    /// Return true if given x values are close enough to each other
148    fn is_converged(&mut self, x1: F, x2: F) -> bool {
149        (x1 - x2).abs() < self.abs()
150    }
151    /// Return true if no more iterations desired
152    fn is_iteration_limit_reached(&mut self, iter: usize) -> bool {
153        iter >= 30
154    }
155}
156
157pub mod brent;
158pub mod eigen;
159pub mod inverse_quadratic;
160pub mod newton_raphson;
161pub mod polynom;
162pub mod regula_falsi;
163pub mod secant;
164
165pub mod debug_convergency;
166pub mod simple_convergency;
167
168#[cfg(test)]
169mod test {
170    use super::*;
171
172    #[test]
173    fn sample_bracketed() {
174        let sample1 = Sample { x: 0f64, y: 0f64 };
175        let sample2 = Sample { x: 1f64, y: 1f64 };
176        let sample3 = Sample { x: 1f64, y: -1f64 };
177        let sample4 = Sample { x: -1f64, y: 0f64 };
178        let sample5 = Sample { x: -1f64, y: 1f64 };
179        assert_eq!(true, sample1.is_bracketed_with(&sample2));
180        assert_eq!(true, sample1.is_bracketed_with(&sample3));
181        assert_eq!(true, sample1.is_bracketed_with(&sample4));
182        assert_eq!(true, sample1.is_bracketed_with(&sample5));
183        assert_eq!(true, sample2.is_bracketed_with(&sample3));
184        assert_eq!(true, sample2.is_bracketed_with(&sample4));
185        assert_eq!(false, sample2.is_bracketed_with(&sample5));
186        assert_eq!(true, sample3.is_bracketed_with(&sample4));
187        assert_eq!(true, sample3.is_bracketed_with(&sample5));
188        assert_eq!(true, sample4.is_bracketed_with(&sample5));
189    }
190
191    #[test]
192    fn root_interval_bracketed() {
193        let sut1 = Interval {
194            begin: Sample { x: 0f64, y: 0f64 },
195            end: Sample { x: 0f64, y: 0f64 },
196        };
197        let sut2 = Interval {
198            begin: Sample { x: 0f32, y: 0f32 },
199            end: Sample { x: 1f32, y: 0f32 },
200        };
201        let sut3 = Interval {
202            begin: Sample { x: 0f64, y: 0f64 },
203            end: Sample { x: 0f64, y: 1f64 },
204        };
205        let sut4 = Interval {
206            begin: Sample { x: -1f64, y: 0f64 },
207            end: Sample { x: 0f64, y: 0f64 },
208        };
209        let sut5 = Interval {
210            begin: Sample { x: -1f64, y: 0f64 },
211            end: Sample { x: 0f64, y: 1f64 },
212        };
213        let sut6 = Interval {
214            begin: Sample { x: -1f32, y: -1f32 },
215            end: Sample { x: 0f32, y: 1f32 },
216        };
217        let sut7 = Interval {
218            begin: Sample { x: 0f64, y: 1f64 },
219            end: Sample { x: 1f64, y: -1f64 },
220        };
221        assert_eq!(true, sut1.is_bracketed());
222        assert_eq!(true, sut2.is_bracketed());
223        assert_eq!(true, sut3.is_bracketed());
224        assert_eq!(true, sut4.is_bracketed());
225        assert_eq!(true, sut5.is_bracketed());
226        assert_eq!(true, sut6.is_bracketed());
227        assert_eq!(true, sut7.is_bracketed());
228    }
229
230    #[test]
231    fn root_interval_not_bracketed() {
232        let sut1 = Interval {
233            begin: Sample { x: 0f64, y: 1f64 },
234            end: Sample { x: 1f64, y: 1f64 },
235        };
236        let sut2 = Interval {
237            begin: Sample { x: -1f64, y: -1f64 },
238            end: Sample { x: 1f64, y: -1f64 },
239        };
240        assert_eq!(false, sut1.is_bracketed());
241        assert_eq!(false, sut2.is_bracketed());
242    }
243
244    #[test]
245    fn root_interval_middle() {
246        let sut1 = Interval {
247            begin: Sample { x: 0f64, y: 1f64 },
248            end: Sample { x: 2f64, y: -3f64 },
249        };
250        let sut2 = Interval {
251            begin: Sample { x: -1f64, y: 0f64 },
252            end: Sample { x: 1f64, y: 0f64 },
253        };
254        assert_eq!(0.5f64, sut1.middle());
255        assert_eq!(0f64, sut2.middle());
256    }
257}
roots/numerical/mod.rs

roots/numerical/
mod.rs
