1 /*****************************************************************************
2 * Copyright (C) 2004-2018 The pagmo development team, *
3 * Advanced Concepts Team (ACT), European Space Agency (ESA) *
4 * http://apps.sourceforge.net/mediawiki/pagmo *
5 * http://apps.sourceforge.net/mediawiki/pagmo/index.php?title=Developers *
6 * http://apps.sourceforge.net/mediawiki/pagmo/index.php?title=Credits *
7 * act@esa.int *
8 * *
9 * This program is free software; you can redistribute it and/or modify *
10 * it under the terms of the GNU General Public License as published by *
11 * the Free Software Foundation; either version 2 of the License, or *
12 * (at your option) any later version. *
13 * *
14 * This program is distributed in the hope that it will be useful, *
15 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
17 * GNU General Public License for more details. *
18 * *
19 * You should have received a copy of the GNU General Public License *
20 * along with this program; if not, write to the *
21 * Free Software Foundation, Inc., *
22 * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
23 *****************************************************************************/
24
25 #ifndef NEWTON_RAPHSON_H
26 #define NEWTON_RAPHSON_H
27 #include <math.h>
28
29 namespace kep_toolbox
30 {
31 /// Newton-Raphson method
32 /**
33 * Standard implementation of the Newton-Raphson method to solve a non-linear equation
34 *
35 * \param[in] x Starting point
36 * \param[in] F equation to be solved in the form F = 0. Needs to be callable as F(x)
37 * \param[in] dF derivative of F with respect to x. Needs to be callable as dF(x)
38 *
39 * @author Dario Izzo (dario.izzo@esa.int)
40 */
41 template <class my_float, class my_functionA, class my_functionB>
newton_raphson(my_float & x,my_functionA F,my_functionB dF,int max_loop,const double & accuracy)42 int newton_raphson(my_float &x, my_functionA F, my_functionB dF, int max_loop, const double &accuracy)
43 {
44 my_float term;
45 // double start = x; //DS: unused
46 // main iteration
47 do {
48 term = (F)(x) / (dF)(x);
49 x = x - term;
50 }
51 // check if term is within required accuracy or loop limit is exceeded
52 while ((fabs(term / std::max(std::fabs(x), 1.)) > accuracy) && (--max_loop));
53 return max_loop;
54 }
55 }
56 /*
57 //----------------------------------------------------------------------------//
58 // test functions
59
60 double func_1(double x) // root is 1.85792
61 { return (cosh(x) + cos(x) - 3.0); } // f(x) = cosh(x) + cos(x) - 3 = 0
62
63 double fdiv_1(double x)
64 { return (sinh(x) - sin(x)); } // f'(x) = sinh(x) - sin(x)
65
66 double func_2(double x) // root is 5.0
67 { return (x*x - 25.0); } // f(x) = x * x - 25 = 0
68
69 double fdiv_2(double x)
70 { return (2.0 * x); } // f'(x) = 2x
71
72 complex<double> func_3(complex<double> x) // roots 5 + or - 3
73 { return x*x - 10.0*x + 34.0; } // f(x) x^2 - 10x + 34
74
75 complex<double> fdiv_3(complex<double> x)
76 { return 2.0*x -10.0; } // f'(x) 2x - 10
77
78
79 struct f3 {
80 complex<double> operator()(complex<double> x) const
81 {
82 return x*x - 10.0*x + 34.0;
83 }
84 };
85
86 struct fd3 {
87 complex<double> operator()(complex<double> x) const
88 {
89 return 2.0*x -10.0;
90 }
91 };
92
93
94 double func_4(double x) // three real roots 4, -3, 1
95 { return 2*x*x*x - 4*x*x - 22*x + 24 ; } // f(x) = 2x^3 - 4x^2 - 22x + 24
96
97 double fdiv_4(double x)
98 { return 6*x*x - 8*x - 22; } // f'(x) = 6x^2 - 8x - 22
99
100 double func_5(double E, double M, double e)
101 { return (E - e * sin(E) - M); }
102
103 double fdiv_5(double E, double e)
104 { return (1 - e * cos(E)); }
105
106 //----------------------------------------------------------------------------//
107 // Main program to test above function
108 int main()
109 {
110 cout << "\nFind root of f(x) = cosh(x) + cos(x) - 3 = 0";
111 double x = 1.0; // initial 'guess' at root
112 if (newton(x, func_1, fdiv_1, 100, 1.0e-8))
113 cout << "\n root x = " << x << ", test of f(x) = " << func_1(x);
114 else cout << "\n failed to find root ";
115
116 cout << "\n\nFind root of f(x) = x * x - 25 = 0";
117 x = 1.0; // initial 'guess' at root
118 if ( newton(x, func_2, fdiv_2, 100, 1.0e-8))
119 cout << "\n root x = " << x << ", test of f(x) = " << func_2(x);
120 else cout << "\n failed to find root ";
121
122 cout << "\n\nFind root of f(x) = x^2 - 10x + 34 = 0";
123 complex<double> xc = complex<double>(1.0, 1.0); // initial 'guess' at root
124 if ( newton(xc, func_3, fdiv_3, 100, 1.0e-8))
125 cout << "\n root x = " << xc << ", test of f(x) = " << func_3(xc);
126 else cout << "\n failed to find root ";
127
128 cout << "\n\nFind root of f(x) = x^2 - 10x + 34 = 0";
129 xc = complex<double>(1.0, -1.0); // initial 'guess' at root
130 if ( newton(xc, func_3, fdiv_3, 100, 1.0e-8))
131 cout << "\n root x = " << xc << ", test of f(x) = " << func_3(xc);
132 else cout << "\n failed to find root ";
133
134 cout << "\n\nFind root of f(x) = 2x^3 - 4x^2 - 22x + 24 = 0";
135 x = 5.0; // initial 'guess' at root
136 if ( newton(x, func_4, fdiv_4, 100, 1.0e-8) )
137 cout << "\n root x = " << x << ", test of f(x) = " << func_4(x);
138 else cout << "\n failed to find root ";
139
140 cout << "\n\nFind root of f(x) = E - e sin(E) - M";
141 x = 0; // initial 'guess' at root
142 newton(x, boost::bind(func_5, _1, 0.2, 0.4), boost::bind(fdiv_5, _1, 0.4), 100, 1.0e-8);
143 if ( newton(x, boost::bind(func_5, _1, 0.2, 0.4), boost::bind(fdiv_5, _1, 0.4), 100, 1.0e-8) )
144 cout << "\n root x = " << x << ", test of f(x) = " << boost::bind(func_5, _1, 0.2,0.4)(x);
145 else cout << "\n failed to find root ";
146 cout << endl;
147 return 0;
148 }
149 */
150
151 #endif // NEWTON_RAPHSON_H
152