KADATH
zerosec.cpp
1 /*
2  Copyright 2017 Philippe Grandclement
3 
4  This file is part of Kadath.
5 
6  Kadath is free software: you can redistribute it and/or modify
7  it under the terms of the GNU General Public License as published by
8  the Free Software Foundation, either version 3 of the License, or
9  (at your option) any later version.
10 
11  Kadath is distributed in the hope that it will be useful,
12  but WITHOUT ANY WARRANTY; without even the implied warranty of
13  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  GNU General Public License for more details.
15 
16  You should have received a copy of the GNU General Public License
17  along with Kadath. If not, see <http://www.gnu.org/licenses/>.
18 */
19 
20 /*
21  * Search for a zero of a function in a given interval, by means of a
22  * secant method.
23  *
24  */
25 
26 // Headers C
27 #include <math.h>
28 #include <assert.h>
29 
30 // Headers Lorene
31 #include "headcpp.hpp"
32 #include "param.hpp"
33 //****************************************************************************
34 namespace Kadath {
35 double zerosec(double (*f)(double, const Param&), const Param& parf,
36  double x1, double x2, double precis, int nitermax, int& niter) {
37 
38  double f0_prec, f0, x0, x0_prec, dx, df ;
39 
40 // Teste si un zero unique existe dans l'intervalle [x_1,x_2]
41 
42 // bool warning = false ;
43 
44  f0_prec = f(x1, parf) ;
45  f0 = f(x2, parf) ;
46  if ( f0*f0_prec > 0.) {
47 // warning = true ;
48  //cout <<
49  // "WARNING: zerosec: there does not exist a unique zero of the function"
50  //<< endl ;
51  //cout << " between x1 = " << x1 << " ( f(x1)=" << f0_prec << " )" << endl ;
52  //cout << " and x2 = " << x2 << " ( f(x2)=" << f0 << " )" << endl ;
53  }
54 
55 // Choisit la borne avec la plus petite valeur de |f(x)| comme la valeur la
56 // "plus recente" de x0
57 
58  if ( fabs(f0) < fabs(f0_prec) ) { // On a bien choisi f0_prec et f0
59  x0_prec = x1 ;
60  x0 = x2 ;
61  }
62  else { // il faut interchanger f0_prec et f0
63  x0_prec = x2 ;
64  x0 = x1 ;
65  double swap = f0_prec ;
66  f0_prec = f0 ;
67  f0 = swap ;
68  }
69 
70 // Debut des iterations de la methode de la secante
71 
72  niter = 0 ;
73  do {
74  df = f0 - f0_prec ;
75  assert(df != double(0)) ;
76  dx = (x0_prec - x0) * f0 / df ;
77  x0_prec = x0 ;
78  f0_prec = f0 ;
79  x0 += dx ;
80  f0 = f(x0, parf) ;
81  niter++ ;
82  if (niter > nitermax) {
83  cout << "zerosec: Maximum number of iterations has been reached ! "
84  << endl ;
85  abort () ;
86  }
87  }
88  while ( ( fabs(dx) > precis ) && ( fabs(f0) > PRECISION ) ) ;
89 
90  // if (warning) {
91 // cout << " A zero has been found at x0 = " << x0 << endl ;
92  // }
93 
94  return x0 ;
95 }
96 
97 
98 }