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
31char 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//****************************************************************************
66namespace Lorene {
67
68double 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}
Time evolution with partial storage (*** under development ***).
Definition evolution.h:371
Parameter storage.
Definition param.h:125
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.
Lorene prototypes.
Definition app_hor.h:64