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::Evolution_std
Time evolution with partial storage (*** under development ***).
Definition
evolution.h:371
Lorene::Param
Parameter storage.
Definition
param.h:125
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
Lorene prototypes.
Definition
app_hor.h:64
C++
Source
Non_class_members
Utilities
zerosec.C
Generated by
1.9.8