LORENE
zerosec_borne.C
1 /*
2  * Search for a zero of a function in a given interval, by means of a
3  * secant method. The input parameters x1 and x2 must be such that
4  * f(x1)*f(x2) < 0 . The obtained root is then necessarily in the
5  * interval [x1,x2].
6  *
7  */
8 
9 /*
10  * Copyright (c) 2002 Jerome Novak
11  *
12  * This file is part of LORENE.
13  *
14  * LORENE is free software; you can redistribute it and/or modify
15  * it under the terms of the GNU General Public License as published by
16  * the Free Software Foundation; either version 2 of the License, or
17  * (at your option) any later version.
18  *
19  * LORENE is distributed in the hope that it will be useful,
20  * but WITHOUT ANY WARRANTY; without even the implied warranty of
21  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
22  * GNU General Public License for more details.
23  *
24  * You should have received a copy of the GNU General Public License
25  * along with LORENE; if not, write to the Free Software
26  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
27  *
28  */
29 
30 
31 
32 
33 /*
34  * $Id: zerosec_borne.C,v 1.8 2016/12/05 16:18:11 j_novak Exp $
35  * $Log: zerosec_borne.C,v $
36  * Revision 1.8 2016/12/05 16:18:11 j_novak
37  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
38  *
39  * Revision 1.7 2014/10/13 08:53:32 j_novak
40  * Lorene classes and functions now belong to the namespace Lorene.
41  *
42  * Revision 1.6 2014/10/06 15:16:11 j_novak
43  * Modified #include directives to use c++ syntax.
44  *
45  * Revision 1.5 2002/10/16 14:37:12 j_novak
46  * Reorganization of #include instructions of standard C++, in order to
47  * use experimental version 3 of gcc.
48  *
49  * Revision 1.4 2002/05/02 15:16:22 j_novak
50  * Added functions for more general bi-fluid EOS
51  *
52  * Revision 1.3 2002/04/18 09:17:17 j_novak
53  * *** empty log message ***
54  *
55  *
56  * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/Utilities/zerosec_borne.C,v 1.8 2016/12/05 16:18:11 j_novak Exp $
57  *
58  */
59 
60 // Headers C
61 #include <cstdlib>
62 #include <cmath>
63 #include <cassert>
64 
65 // Headers Lorene
66 #include "headcpp.h"
67 #include "param.h"
68 //****************************************************************************
69 namespace Lorene {
70 
71 double zerosec_b(double (*f)(double, const Param&), const Param& parf,
72  double x1, double x2, double precis, int nitermax, int& niter) {
73 
74  double f0_moins, f0, x0, x0_moins, dx, df , fnew, xnew;
75 
76 // Teste si un zero unique existe dans l'intervalle [x_1,x_2]
77 
78  f0_moins = f(x1, parf) ;
79  f0 = f(x2, parf) ;
80  if ( f0*f0_moins > 0.) {
81  cout <<
82  "WARNING: zerosec: there may not exist a zero of the function"
83  << endl ;
84  cout << " between x1 = " << x1 << " ( f(x1)=" << f0_moins << " )" << endl ;
85  cout << " and x2 = " << x2 << " ( f(x2)=" << f0 << " )" << endl ;
86  }
87 
88 // Choisit la borne avec la plus grande valeur de f(x) (borne positive)
89 // comme la valeur la de x0
90 
91  if ( f0_moins < f0) { // On a bien choisi f0_moins et f0
92  x0_moins = x1 ;
93  x0 = x2 ;
94  }
95  else { // il faut interchanger f0_moins et f0
96  x0_moins = x2 ;
97  x0 = x1 ;
98  double swap = f0_moins ;
99  f0_moins = f0 ;
100  f0 = swap ;
101  }
102 
103 // Debut des iterations de la methode de la secante
104 
105  niter = 0 ;
106  do {
107  df = f0 - f0_moins ;
108  assert(df != double(0)) ;
109  xnew = (x0_moins*f0 - x0*f0_moins)/df ; ;
110  fnew = f(xnew, parf) ;
111  if (fnew < 0.) {
112  dx = x0_moins - xnew ;
113  x0_moins = xnew ;
114  f0_moins = fnew ;
115  }
116  else {
117  dx = x0 - xnew ;
118  x0 = xnew ;
119  f0 = fnew ;
120  }
121  niter++ ;
122  if (niter > nitermax) {
123 //cout << "zerosec: Maximum number of iterations has been reached ! "
124  // << endl ;
125  //cout << x0_moins << ", " << xnew << ", " << x0 << endl ;
126  //cout << f0_moins << ", " << fnew << ", " << f0 << endl ;
127 
128  return xnew ;
129  }
130  }
131  while ( ( fabs(dx) > precis ) && ( fabs(fnew) > 1.e-15 ) ) ;
132 
133  return xnew ;
134 }
135 
136 
137 
138 }
Lorene prototypes.
Definition: app_hor.h:67
double zerosec_b(double(*f)(double, const Param &), const Param &par, double a, double b, double precis, int nitermax, int &niter)
Finding the zero a function on a bounded domain.
Definition: zerosec_borne.C:71
Parameter storage.
Definition: param.h:125