1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
|
/*
Copyright (C) 2000, 2001, 2002 RiskMap srl
This file is part of QuantLib, a free-software/open-source library
for financial quantitative analysts and developers - http://quantlib.org/
QuantLib is free software: you can redistribute it and/or modify it under the
terms of the QuantLib license. You should have received a copy of the
license along with this program; if not, please email ferdinando@ametrano.net
The license is also available online at http://quantlib.org/html/license.html
This program is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
FOR A PARTICULAR PURPOSE. See the license for more details.
*/
/*! \file brent.cpp
\brief Brent 1-D solver
\fullpath
ql/Solvers1D/%brent.cpp
*/
// $Id: brent.cpp,v 1.6 2002/01/16 14:41:17 nando Exp $
/* The implementation of the algorithm was inspired by
* "Numerical Recipes in C", 2nd edition, Press, Teukolsky, Vetterling, Flannery
* Chapter 9
*/
#include <ql/Solvers1D/brent.hpp>
namespace QuantLib {
namespace Solvers1D {
#define SIGN(a,b) ((b) >= 0.0 ? QL_FABS(a) : -QL_FABS(a))
double Brent::solve_(const ObjectiveFunction& f,
double xAccuracy) const {
double min1, min2;
double froot, p, q, r, s, xAcc1, xMid;
// dummy assignements to avoid compiler warning
double d = 0.0;
double e = 0.0;
root_ = xMax_;
froot = fxMax_;
while (evaluationNumber_<=maxEvaluations_) {
if ((froot > 0.0 && fxMax_ > 0.0) ||
(froot < 0.0 && fxMax_ < 0.0)) {
// Rename xMin_, root_, xMax_ and adjust bounding interval d
xMax_=xMin_;
fxMax_=fxMin_;
e=d=root_-xMin_;
}
if (QL_FABS(fxMax_) < QL_FABS(froot)) {
xMin_=root_;
root_=xMax_;
xMax_=xMin_;
fxMin_=froot;
froot=fxMax_;
fxMax_=fxMin_;
}
// Convergence check
xAcc1=2.0*QL_EPSILON*QL_FABS(root_)+0.5*xAccuracy;
xMid=(xMax_-root_)/2.0;
if (QL_FABS(xMid) <= xAcc1 || froot == 0.0) return root_;
if (QL_FABS(e) >= xAcc1 && QL_FABS(fxMin_) > QL_FABS(froot)) {
s=froot/fxMin_; // Attempt inverse quadratic interpolation
if (xMin_ == xMax_) {
p=2.0*xMid*s;
q=1.0-s;
} else {
q=fxMin_/fxMax_;
r=froot/fxMax_;
p=s*(2.0*xMid*q*(q-r)-(root_-xMin_)*(r-1.0));
q=(q-1.0)*(r-1.0)*(s-1.0);
}
if (p > 0.0) q = -q; // Check whether in bounds
p=QL_FABS(p);
min1=3.0*xMid*q-QL_FABS(xAcc1*q);
min2=QL_FABS(e*q);
if (2.0*p < (min1 < min2 ? min1 : min2)) {
e=d; // Accept interpolation
d=p/q;
} else {
d=xMid; // Interpolation failed, use bisection
e=d;
}
} else { // Bounds decreasing too slowly, use bisection
d=xMid;
e=d;
}
xMin_=root_;
fxMin_=froot;
if (QL_FABS(d) > xAcc1)
root_ += d;
else
root_ += SIGN(xAcc1,xMid);
froot=f(root_);
evaluationNumber_++;
}
throw Error("Brent: maximum number of function evaluations ("
+ IntegerFormatter::toString(maxEvaluations_) + ") exceeded");
QL_DUMMY_RETURN(0.0);
}
}
}
|
