LORENE
zerosec.C
1 /*
2  * Search for a zero of a function in a given interval, by means of a
3  * secant method.
4  *
5  */
6 
7 /*
8  * Copyright (c) 1999-2001 Eric Gourgoulhon
9  *
10  * This file is part of LORENE.
11  *
12  * LORENE is free software; you can redistribute it and/or modify
13  * it under the terms of the GNU General Public License as published by
14  * the Free Software Foundation; either version 2 of the License, or
15  * (at your option) any later version.
16  *
17  * LORENE is distributed in the hope that it will be useful,
18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20  * GNU General Public License for more details.
21  *
22  * You should have received a copy of the GNU General Public License
23  * along with LORENE; if not, write to the Free Software
24  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
25  *
26  */
27 
28 
29 
30 
31 /*
32  * $Id: zerosec.C,v 1.7 2016/12/05 16:18:11 j_novak Exp $
33  * $Log: zerosec.C,v $
34  * Revision 1.7 2016/12/05 16:18:11 j_novak
35  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
36  *
37  * Revision 1.6 2014/10/13 08:53:32 j_novak
38  * Lorene classes and functions now belong to the namespace Lorene.
39  *
40  * Revision 1.5 2014/07/04 12:09:06 j_novak
41  * New argument in zerosec(): a boolean (false by default) for aborting if the number of iteration is greater than the max.
42  *
43  * Revision 1.4 2002/10/16 14:37:12 j_novak
44  * Reorganization of #include instructions of standard C++, in order to
45  * use experimental version 3 of gcc.
46  *
47  * Revision 1.3 2002/04/11 09:19:46 j_novak
48  * Back to old version of zerosec
49  *
50  * Revision 1.1.1.1 2001/11/20 15:19:29 e_gourgoulhon
51  * LORENE
52  *
53  * Revision 1.6 2001/10/17 08:16:47 eric
54  * In case there is not a single zero in the interval, the found
55  * zero is displayed in the warning message.
56  *
57  * Revision 1.5 2000/01/04 13:20:34 eric
58  * Test final f0 != double(0) remplace par fabs(f0) > 1.e-15 .
59  *
60  * Revision 1.4 1999/12/20 09:46:08 eric
61  * Anglicisation des messages.
62  *
63  * Revision 1.3 1999/12/17 10:08:46 eric
64  * Le test final fabs(f0) > 1.e-14 est remplace par f0 != 0.
65  *
66  * Revision 1.2 1999/12/17 09:37:40 eric
67  * Ajout de assert(df != 0).
68  *
69  * Revision 1.1 1999/12/15 09:41:34 eric
70  * Initial revision
71  *
72  *
73  * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/Utilities/zerosec.C,v 1.7 2016/12/05 16:18:11 j_novak Exp $
74  *
75  */
76 
77 // Headers C
78 #include <cstdlib>
79 #include <cmath>
80 #include <cassert>
81 
82 // Headers C++
83 #include <exception>
84 
85 // Headers Lorene
86 #include "headcpp.h"
87 #include "param.h"
88 //****************************************************************************
89 
90 namespace Lorene {
91 
92 double zerosec(double (*f)(double, const Param&), const Param& parf,
93  double x1, double x2, double precis, int nitermax, int& niter,
94  bool abor) {
95 
96  double f0_prec, f0, x0, x0_prec, dx, df ;
97 
98 // Teste si un zero unique existe dans l'intervalle [x_1,x_2]
99 
100  bool warning = false ;
101 
102  f0_prec = f(x1, parf) ;
103  f0 = f(x2, parf) ;
104  if ( f0*f0_prec > 0.) {
105  warning = true ;
106  cout <<
107  "WARNING: zerosec: there does not exist a unique zero of the function"
108  << endl ;
109  cout << " between x1 = " << x1 << " ( f(x1)=" << f0_prec << " )" << endl ;
110  cout << " and x2 = " << x2 << " ( f(x2)=" << f0 << " )" << endl ;
111  }
112 
113 // Choisit la borne avec la plus petite valeur de |f(x)| comme la valeur la
114 // "plus recente" de x0
115 
116  if ( fabs(f0) < fabs(f0_prec) ) { // On a bien choisi f0_prec et f0
117  x0_prec = x1 ;
118  x0 = x2 ;
119  }
120  else { // il faut interchanger f0_prec et f0
121  x0_prec = x2 ;
122  x0 = x1 ;
123  double swap = f0_prec ;
124  f0_prec = f0 ;
125  f0 = swap ;
126  }
127 
128 // Debut des iterations de la methode de la secante
129 
130  niter = 0 ;
131  do {
132  df = f0 - f0_prec ;
133  assert(df != double(0)) ;
134  dx = (x0_prec - x0) * f0 / df ;
135  x0_prec = x0 ;
136  f0_prec = f0 ;
137  x0 += dx ;
138  f0 = f(x0, parf) ;
139  niter++ ;
140  if (niter > nitermax) {
141  cout << "zerosec: Maximum number of iterations has been reached ! "
142  << endl ;
143  if (abor)
144  abort () ;
145  else {
146  warning = true ;
147  f0 = 0. ;
148  }
149  }
150  }
151  while ( ( fabs(dx) > precis ) && ( fabs(f0) > 1.e-15 ) ) ;
152 
153  if (warning) {
154  cout << " A zero may have been found at x0 = " << x0 << endl ;
155  }
156 
157  return x0 ;
158 }
159 
160 
161 
162 }
Lorene prototypes.
Definition: app_hor.h:67
double zerosec(double(*f)(double, const Param &), const Param &par, double a, double b, double precis, int nitermax, int &niter, bool abort=true)
Finding the zero a function.
Definition: zerosec.C:92
Parameter storage.
Definition: param.h:125