KADATH
zerosec.cpp
1
/*
2
Copyright 2017 Philippe Grandclement
3
4
This file is part of Kadath.
5
6
Kadath is free software: you can redistribute it and/or modify
7
it under the terms of the GNU General Public License as published by
8
the Free Software Foundation, either version 3 of the License, or
9
(at your option) any later version.
10
11
Kadath is distributed in the hope that it will be useful,
12
but WITHOUT ANY WARRANTY; without even the implied warranty of
13
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14
GNU General Public License for more details.
15
16
You should have received a copy of the GNU General Public License
17
along with Kadath. If not, see <http://www.gnu.org/licenses/>.
18
*/
19
20
/*
21
* Search for a zero of a function in a given interval, by means of a
22
* secant method.
23
*
24
*/
25
26
// Headers C
27
#include <math.h>
28
#include <assert.h>
29
30
// Headers Lorene
31
#include "headcpp.hpp"
32
#include "param.hpp"
33
//****************************************************************************
34
namespace
Kadath {
35
double
zerosec(
double
(*f)(
double
,
const
Param&),
const
Param& parf,
36
double
x1,
double
x2,
double
precis,
int
nitermax,
int
& niter) {
37
38
double
f0_prec, f0, x0, x0_prec, dx, df ;
39
40
// Teste si un zero unique existe dans l'intervalle [x_1,x_2]
41
42
// bool warning = false ;
43
44
f0_prec = f(x1, parf) ;
45
f0 = f(x2, parf) ;
46
if
( f0*f0_prec > 0.) {
47
// warning = true ;
48
//cout <<
49
// "WARNING: zerosec: there does not exist a unique zero of the function"
50
//<< endl ;
51
//cout << " between x1 = " << x1 << " ( f(x1)=" << f0_prec << " )" << endl ;
52
//cout << " and x2 = " << x2 << " ( f(x2)=" << f0 << " )" << endl ;
53
}
54
55
// Choisit la borne avec la plus petite valeur de |f(x)| comme la valeur la
56
// "plus recente" de x0
57
58
if
( fabs(f0) < fabs(f0_prec) ) {
// On a bien choisi f0_prec et f0
59
x0_prec = x1 ;
60
x0 = x2 ;
61
}
62
else
{
// il faut interchanger f0_prec et f0
63
x0_prec = x2 ;
64
x0 = x1 ;
65
double
swap = f0_prec ;
66
f0_prec = f0 ;
67
f0 = swap ;
68
}
69
70
// Debut des iterations de la methode de la secante
71
72
niter = 0 ;
73
do
{
74
df = f0 - f0_prec ;
75
assert(df !=
double
(0)) ;
76
dx = (x0_prec - x0) * f0 / df ;
77
x0_prec = x0 ;
78
f0_prec = f0 ;
79
x0 += dx ;
80
f0 = f(x0, parf) ;
81
niter++ ;
82
if
(niter > nitermax) {
83
cout <<
"zerosec: Maximum number of iterations has been reached ! "
84
<< endl ;
85
abort () ;
86
}
87
}
88
while
( ( fabs(dx) > precis ) && ( fabs(f0) > PRECISION ) ) ;
89
90
// if (warning) {
91
// cout << " A zero has been found at x0 = " << x0 << endl ;
92
// }
93
94
return
x0 ;
95
}
96
97
98
}
src
Utilities
zerosec.cpp
Generated by
1.9.1