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
C++
Source
Non_class_members
Utilities
zerosec.C
Generated by
1.8.17