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 char zerosec_borne_C[] = "$Header: /cvsroot/Lorene/C++/Source/Non_class_members/Utilities/zerosec_borne.C,v 1.7 2014/10/13 08:53:32 j_novak Exp $" ;
32 
33 /*
34  * $Id: zerosec_borne.C,v 1.7 2014/10/13 08:53:32 j_novak Exp $
35  * $Log: zerosec_borne.C,v $
36  * Revision 1.7 2014/10/13 08:53:32 j_novak
37  * Lorene classes and functions now belong to the namespace Lorene.
38  *
39  * Revision 1.6 2014/10/06 15:16:11 j_novak
40  * Modified #include directives to use c++ syntax.
41  *
42  * Revision 1.5 2002/10/16 14:37:12 j_novak
43  * Reorganization of #include instructions of standard C++, in order to
44  * use experimental version 3 of gcc.
45  *
46  * Revision 1.4 2002/05/02 15:16:22 j_novak
47  * Added functions for more general bi-fluid EOS
48  *
49  * Revision 1.3 2002/04/18 09:17:17 j_novak
50  * *** empty log message ***
51  *
52  *
53  * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/Utilities/zerosec_borne.C,v 1.7 2014/10/13 08:53:32 j_novak Exp $
54  *
55  */
56 
57 // Headers C
58 #include <cstdlib>
59 #include <cmath>
60 #include <cassert>
61 
62 // Headers Lorene
63 #include "headcpp.h"
64 #include "param.h"
65 //****************************************************************************
66 namespace Lorene {
67 
68 double zerosec_b(double (*f)(double, const Param&), const Param& parf,
69  double x1, double x2, double precis, int nitermax, int& niter) {
70 
71  double f0_moins, f0, x0, x0_moins, dx, df , fnew, xnew;
72 
73 // Teste si un zero unique existe dans l'intervalle [x_1,x_2]
74 
75  f0_moins = f(x1, parf) ;
76  f0 = f(x2, parf) ;
77  if ( f0*f0_moins > 0.) {
78  cout <<
79  "WARNING: zerosec: there may not exist a zero of the function"
80  << endl ;
81  cout << " between x1 = " << x1 << " ( f(x1)=" << f0_moins << " )" << endl ;
82  cout << " and x2 = " << x2 << " ( f(x2)=" << f0 << " )" << endl ;
83  }
84 
85 // Choisit la borne avec la plus grande valeur de f(x) (borne positive)
86 // comme la valeur la de x0
87 
88  if ( f0_moins < f0) { // On a bien choisi f0_moins et f0
89  x0_moins = x1 ;
90  x0 = x2 ;
91  }
92  else { // il faut interchanger f0_moins et f0
93  x0_moins = x2 ;
94  x0 = x1 ;
95  double swap = f0_moins ;
96  f0_moins = f0 ;
97  f0 = swap ;
98  }
99 
100 // Debut des iterations de la methode de la secante
101 
102  niter = 0 ;
103  do {
104  df = f0 - f0_moins ;
105  assert(df != double(0)) ;
106  xnew = (x0_moins*f0 - x0*f0_moins)/df ; ;
107  fnew = f(xnew, parf) ;
108  if (fnew < 0.) {
109  dx = x0_moins - xnew ;
110  x0_moins = xnew ;
111  f0_moins = fnew ;
112  }
113  else {
114  dx = x0 - xnew ;
115  x0 = xnew ;
116  f0 = fnew ;
117  }
118  niter++ ;
119  if (niter > nitermax) {
120 //cout << "zerosec: Maximum number of iterations has been reached ! "
121  // << endl ;
122  //cout << x0_moins << ", " << xnew << ", " << x0 << endl ;
123  //cout << f0_moins << ", " << fnew << ", " << f0 << endl ;
124 
125  return xnew ;
126  }
127  }
128  while ( ( fabs(dx) > precis ) && ( fabs(fnew) > 1.e-15 ) ) ;
129 
130  return xnew ;
131 }
132 
133 
134 
135 }
Lorene
Lorene prototypes.
Definition: app_hor.h:64
Lorene::zerosec_b
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:68
Lorene::Param
Parameter storage.
Definition: param.h:125