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 char zerosec_C[] = "$Header: /cvsroot/Lorene/C++/Source/Non_class_members/Utilities/zerosec.C,v 1.6 2014/10/13 08:53:32 j_novak Exp $" ;
30 
31 /*
32  * $Id: zerosec.C,v 1.6 2014/10/13 08:53:32 j_novak Exp $
33  * $Log: zerosec.C,v $
34  * Revision 1.6 2014/10/13 08:53:32 j_novak
35  * Lorene classes and functions now belong to the namespace Lorene.
36  *
37  * Revision 1.5 2014/07/04 12:09:06 j_novak
38  * New argument in zerosec(): a boolean (false by default) for aborting if the number of iteration is greater than the max.
39  *
40  * Revision 1.4 2002/10/16 14:37:12 j_novak
41  * Reorganization of #include instructions of standard C++, in order to
42  * use experimental version 3 of gcc.
43  *
44  * Revision 1.3 2002/04/11 09:19:46 j_novak
45  * Back to old version of zerosec
46  *
47  * Revision 1.1.1.1 2001/11/20 15:19:29 e_gourgoulhon
48  * LORENE
49  *
50  * Revision 1.6 2001/10/17 08:16:47 eric
51  * In case there is not a single zero in the interval, the found
52  * zero is displayed in the warning message.
53  *
54  * Revision 1.5 2000/01/04 13:20:34 eric
55  * Test final f0 != double(0) remplace par fabs(f0) > 1.e-15 .
56  *
57  * Revision 1.4 1999/12/20 09:46:08 eric
58  * Anglicisation des messages.
59  *
60  * Revision 1.3 1999/12/17 10:08:46 eric
61  * Le test final fabs(f0) > 1.e-14 est remplace par f0 != 0.
62  *
63  * Revision 1.2 1999/12/17 09:37:40 eric
64  * Ajout de assert(df != 0).
65  *
66  * Revision 1.1 1999/12/15 09:41:34 eric
67  * Initial revision
68  *
69  *
70  * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/Utilities/zerosec.C,v 1.6 2014/10/13 08:53:32 j_novak Exp $
71  *
72  */
73 
74 // Headers C
75 #include <cstdlib>
76 #include <cmath>
77 #include <cassert>
78 
79 // Headers C++
80 #include <exception>
81 
82 // Headers Lorene
83 #include "headcpp.h"
84 #include "param.h"
85 //****************************************************************************
86 
87 namespace Lorene {
88 
89 double zerosec(double (*f)(double, const Param&), const Param& parf,
90  double x1, double x2, double precis, int nitermax, int& niter,
91  bool abor) {
92 
93  double f0_prec, f0, x0, x0_prec, dx, df ;
94 
95 // Teste si un zero unique existe dans l'intervalle [x_1,x_2]
96 
97  bool warning = false ;
98 
99  f0_prec = f(x1, parf) ;
100  f0 = f(x2, parf) ;
101  if ( f0*f0_prec > 0.) {
102  warning = true ;
103  cout <<
104  "WARNING: zerosec: there does not exist a unique zero of the function"
105  << endl ;
106  cout << " between x1 = " << x1 << " ( f(x1)=" << f0_prec << " )" << endl ;
107  cout << " and x2 = " << x2 << " ( f(x2)=" << f0 << " )" << endl ;
108  }
109 
110 // Choisit la borne avec la plus petite valeur de |f(x)| comme la valeur la
111 // "plus recente" de x0
112 
113  if ( fabs(f0) < fabs(f0_prec) ) { // On a bien choisi f0_prec et f0
114  x0_prec = x1 ;
115  x0 = x2 ;
116  }
117  else { // il faut interchanger f0_prec et f0
118  x0_prec = x2 ;
119  x0 = x1 ;
120  double swap = f0_prec ;
121  f0_prec = f0 ;
122  f0 = swap ;
123  }
124 
125 // Debut des iterations de la methode de la secante
126 
127  niter = 0 ;
128  do {
129  df = f0 - f0_prec ;
130  assert(df != double(0)) ;
131  dx = (x0_prec - x0) * f0 / df ;
132  x0_prec = x0 ;
133  f0_prec = f0 ;
134  x0 += dx ;
135  f0 = f(x0, parf) ;
136  niter++ ;
137  if (niter > nitermax) {
138  cout << "zerosec: Maximum number of iterations has been reached ! "
139  << endl ;
140  if (abor)
141  abort () ;
142  else {
143  warning = true ;
144  f0 = 0. ;
145  }
146  }
147  }
148  while ( ( fabs(dx) > precis ) && ( fabs(f0) > 1.e-15 ) ) ;
149 
150  if (warning) {
151  cout << " A zero may have been found at x0 = " << x0 << endl ;
152  }
153 
154  return x0 ;
155 }
156 
157 
158 
159 }
Lorene
Lorene prototypes.
Definition: app_hor.h:64
Lorene::zerosec
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:89
Lorene::Param
Parameter storage.
Definition: param.h:125