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
29char 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
87namespace Lorene {
88
89double 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}
Time evolution with partial storage (*** under development ***).
Definition evolution.h:371
Parameter storage.
Definition param.h:125
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 prototypes.
Definition app_hor.h:64