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
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112

pub use crate::traits::{Number, Signed, Zero, One};
pub use crate::complex::Complex;
pub use crate::vector::{Vector, Vec64};
pub use crate::matrix::{Matrix, Mat64};

pub struct Newton<T> {
    tol: f64,
    delta: f64,
    max_iter: usize,
    guess: T,
}

impl<T> Newton<T> {
    /// Create a new newton object
    #[inline]
    pub const fn new( guess: T ) -> Self {
        let tol: f64 = 1.0e-8;
        let delta: f64 = 1.0e-8;
        let max_iter: usize = 20;
        Newton { tol, delta, max_iter, guess }
    }

    /// Edit the convergence tolerance 
    #[inline]
    pub fn tolerance(&mut self, tolerance: f64 ) {
        self.tol = tolerance;
    }

    /// Edit the finite difference derivative step 
    #[inline]
    pub fn delta(&mut self, delta: f64 ) {
        self.delta = delta;
    }

    /// Edit the maximum number of iterations 
    #[inline]
    pub fn iterations(&mut self, iterations: usize ) {
        self.max_iter = iterations;
    }

    /// Edit the initial guess 
    #[inline]
    pub fn guess(&mut self, guess: T ) {
        self.guess = guess;
    }

}

impl<T: Copy> Newton<T> {
    /// Return the current parameters as a tuple 
    #[inline]
    pub fn parameters(&self) -> ( f64, f64, usize, T ) {
        let parameters = ( self.tol, self.delta, self.max_iter, self.guess );
        parameters
    }
}

impl Newton<f64> {
    /// Solve the equation via Newton iteration 
    #[inline]
    pub fn solve(&self, func: &dyn Fn(f64) -> f64 ) -> Result<f64, f64> {
        let mut current: f64 = self.guess;
        for _ in 0..self.max_iter {
            let deriv = ( func( current + self.delta ) - 
                          func( current - self.delta ) ) / ( 2.0 * self.delta );
            let dx = func(current) / deriv;
            current -= dx;
            if dx.abs() <= self.tol {
                return Ok( current );
            }
        }
        Err( current ) 
    }
}

impl Newton<Vec64> {
    /// Solve the vector equation via Newton iteration 
    #[inline] 
    pub fn solve(&self, func: &dyn Fn(Vec64) -> Vec64) -> Result<Vec64, Vec64> {
        let mut current: Vec64 = self.guess.clone();
        for _ in 0..self.max_iter {
            let f: Vec64 = func( current.clone() );
            let max_residual = f.norm_inf();
            let mut j = Mat64::jacobian( current.clone(), func, self.delta );
            let dx: Vec64 = j.solve_basic( f );
            current -= dx;
            if max_residual <= self.tol {
                return Ok( current )
            }
        }
        Err( current )
    }

    /// Solve the vector equation via Newton iteration using the exact Jacobian
    #[inline] 
    pub fn solve_jacobian(&self, func: &dyn Fn(Vec64) -> Vec64, 
                                 jac: &dyn Fn(Vec64) -> Mat64 ) -> Result<Vec64, Vec64> {
        let mut current: Vec64 = self.guess.clone();
        for _ in 0..self.max_iter {
            let f: Vec64 = func( current.clone() );
            let max_residual = f.norm_inf();
            let mut j: Mat64 = jac( current.clone() ); 
            let dx: Vec64 = j.solve_basic( f );
            current -= dx;
            if max_residual <= self.tol {
                return Ok( current )
            }
        }
        Err( current )
    }

}