1/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2
3/*
4 Copyright (C) 2003 RiskMap srl
5 Copyright (C) 2012 StatPro Italia srl
6
7 This file is part of QuantLib, a free-software/open-source library
8 for financial quantitative analysts and developers - http://quantlib.org/
9
10 QuantLib is free software: you can redistribute it and/or modify it
11 under the terms of the QuantLib license. You should have received a
12 copy of the license along with this program; if not, please email
13 <quantlib-dev@lists.sf.net>. The license is also available online at
14 <http://quantlib.org/license.shtml>.
15
16 This program is distributed in the hope that it will be useful, but WITHOUT
17 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
18 FOR A PARTICULAR PURPOSE. See the license for more details.
19*/
20
21#include "solvers.hpp"
22#include "utilities.hpp"
23#include <ql/math/solvers1d/brent.hpp>
24#include <ql/math/solvers1d/bisection.hpp>
25#include <ql/math/solvers1d/falseposition.hpp>
26#include <ql/math/solvers1d/ridder.hpp>
27#include <ql/math/solvers1d/secant.hpp>
28#include <ql/math/solvers1d/newton.hpp>
29#include <ql/math/solvers1d/newtonsafe.hpp>
30#include <ql/math/solvers1d/finitedifferencenewtonsafe.hpp>
31
32using namespace QuantLib;
33using namespace boost::unit_test_framework;
34
35namespace {
36
37 class F1 {
38 public:
39 Real operator()(Real x) const { return x*x-1.0; }
40 Real derivative(Real x) const { return 2.0*x; }
41 };
42
43 class F2 {
44 public:
45 Real operator()(Real x) const { return 1.0-x*x; }
46 Real derivative(Real x) const { return -2.0*x; }
47 };
48
49 class F3 {
50 public:
51 Real operator()(Real x) const { return std::atan(x: x-1); }
52 Real derivative(Real x) const { return 1.0 / (1.0+(x-1.0)*(x-1.0)); }
53 };
54
55 template <class S, class F>
56 void test_not_bracketed(const S& solver, const std::string& name,
57 const F& f, Real guess) {
58 Real accuracy[] = { 1.0e-4, 1.0e-6, 1.0e-8 };
59 Real expected = 1.0;
60 for (Real& i : accuracy) {
61 Real root = solver.solve(f, i, guess, 0.1);
62 if (std::fabs(x: root - expected) > i) {
63 BOOST_FAIL(name << " solver (not bracketed):\n"
64 << " expected: " << expected << "\n"
65 << " calculated: " << root << "\n"
66 << " accuracy: " << i);
67 }
68 }
69 }
70
71 template <class S, class F>
72 void test_bracketed(const S& solver, const std::string& name,
73 const F& f, Real guess) {
74 Real accuracy[] = { 1.0e-4, 1.0e-6, 1.0e-8 };
75 Real expected = 1.0;
76 for (Real& i : accuracy) {
77 // guess on the left side of the root, increasing function
78 Real root = solver.solve(f, i, guess, 0.0, 2.0);
79 if (std::fabs(x: root - expected) > i) {
80 BOOST_FAIL(name << " solver (bracketed):\n"
81 << " expected: " << expected << "\n"
82 << " calculated: " << root << "\n"
83 << " accuracy: " << i);
84 }
85 }
86 }
87
88 class Probe {
89 public:
90 Probe(Real& result, Real offset)
91 : result_(result), previous_(result), offset_(offset) {}
92 Real operator()(Real x) const {
93 result_ = x;
94 return previous_ + offset_ - x*x;
95 }
96 Real derivative(Real x) const { return 2.0*x; }
97 private:
98 Real& result_;
99 Real previous_;
100 Real offset_;
101 };
102
103 template <class S>
104 void test_last_call_with_root(const S& solver, const std::string& name,
105 bool bracketed, Real accuracy) {
106
107 Real mins[] = { 3.0, 2.25, 1.5, 1.0 };
108 Real maxs[] = { 7.0, 5.75, 4.5, 3.0 };
109 Real steps[] = { 0.2, 0.2, 0.1, 0.1 };
110 Real offsets[] = { 25.0, 11.0, 5.0, 1.0 };
111 Real guesses[] = { 4.5, 4.5, 2.5, 2.5 };
112 //Real expected[] = { 5.0, 4.0, 3.0, 2.0 };
113
114 Real argument = 0.0;
115 Real result;
116
117 for (Size i=0; i<4; ++i) {
118 if (bracketed) {
119 result = solver.solve(Probe(argument, offsets[i]), accuracy,
120 guesses[i], mins[i], maxs[i]);
121 } else {
122 result = solver.solve(Probe(argument, offsets[i]), accuracy,
123 guesses[i], steps[i]);
124 }
125
126 Real error = std::fabs(x: result-argument);
127 // the solver should have called the function with
128 // the very same value it's returning. But the internal
129 // 80bit length of the x87 FPU register might lead to
130 // a very small glitch when compiled with -mfpmath=387 on gcc
131 if (error > 2*QL_EPSILON) {
132 BOOST_FAIL(name << " solver ("
133 << (bracketed ? "" : "not ")
134 << "bracketed):\n"
135 << " index: " << i << "\n"
136 << " expected: " << result << "\n"
137 << " calculated: " << argument << "\n"
138 << " error: " << error);
139 }
140 }
141 }
142
143 template <class S>
144 void test_solver(const S& solver, const std::string& name, Real accuracy) {
145 // guess on the left side of the root, increasing function
146 test_not_bracketed(solver, name, F1(), 0.5);
147 test_bracketed(solver, name, F1(), 0.5);
148 // guess on the right side of the root, increasing function
149 test_not_bracketed(solver, name, F1(), 1.5);
150 test_bracketed(solver, name, F1(), 1.5);
151 // guess on the left side of the root, decreasing function
152 test_not_bracketed(solver, name, F2(), 0.5);
153 test_bracketed(solver, name, F2(), 0.5);
154 // guess on the right side of the root, decreasing function
155 test_not_bracketed(solver, name, F2(), 1.5);
156 test_bracketed(solver, name, F2(), 1.5);
157 // situation where bisection is used in the finite difference
158 // newton solver as the first step and where the initial
159 // guess is equal to the next estimate (which causes an infinite
160 // derivative if we do not handle this case with special care)
161 test_not_bracketed(solver, name, F3(), 1.00001);
162 // check that the last function call is made with the root value
163 if(accuracy != Null<Real>()) {
164 test_last_call_with_root(solver, name, false, accuracy);
165 test_last_call_with_root(solver, name, true, accuracy);
166 }
167 }
168
169}
170
171
172void Solver1DTest::testBrent() {
173 BOOST_TEST_MESSAGE("Testing Brent solver...");
174 test_solver(solver: Brent(), name: "Brent", accuracy: 1.0e-6);
175}
176
177void Solver1DTest::testBisection() {
178 BOOST_TEST_MESSAGE("Testing bisection solver...");
179 test_solver(solver: Bisection(), name: "Bisection", accuracy: 1.0e-6);
180}
181
182void Solver1DTest::testFalsePosition() {
183 BOOST_TEST_MESSAGE("Testing false-position solver...");
184 test_solver(solver: FalsePosition(), name: "FalsePosition", accuracy: 1.0e-6);
185}
186
187void Solver1DTest::testNewton() {
188 BOOST_TEST_MESSAGE("Testing Newton solver...");
189 test_solver(solver: Newton(), name: "Newton", accuracy: 1.0e-12);
190}
191
192void Solver1DTest::testNewtonSafe() {
193 BOOST_TEST_MESSAGE("Testing Newton-safe solver...");
194 test_solver(solver: NewtonSafe(), name: "NewtonSafe", accuracy: 1.0e-9);
195}
196
197void Solver1DTest::testFiniteDifferenceNewtonSafe() {
198 BOOST_TEST_MESSAGE("Testing finite-difference Newton-safe solver...");
199 test_solver(solver: FiniteDifferenceNewtonSafe(), name: "FiniteDifferenceNewtonSafe", accuracy: Null<Real>());
200}
201
202void Solver1DTest::testRidder() {
203 BOOST_TEST_MESSAGE("Testing Ridder solver...");
204 test_solver(solver: Ridder(), name: "Ridder", accuracy: 1.0e-6);
205}
206
207void Solver1DTest::testSecant() {
208 BOOST_TEST_MESSAGE("Testing secant solver...");
209 test_solver(solver: Secant(), name: "Secant", accuracy: 1.0e-6);
210}
211
212
213test_suite* Solver1DTest::suite() {
214 auto* suite = BOOST_TEST_SUITE("1-D solver tests");
215 suite->add(QUANTLIB_TEST_CASE(&Solver1DTest::testBrent));
216 suite->add(QUANTLIB_TEST_CASE(&Solver1DTest::testBisection));
217 suite->add(QUANTLIB_TEST_CASE(&Solver1DTest::testFalsePosition));
218 suite->add(QUANTLIB_TEST_CASE(&Solver1DTest::testNewton));
219 suite->add(QUANTLIB_TEST_CASE(&Solver1DTest::testNewtonSafe));
220 suite->add(QUANTLIB_TEST_CASE(&Solver1DTest::testFiniteDifferenceNewtonSafe));
221 suite->add(QUANTLIB_TEST_CASE(&Solver1DTest::testRidder));
222 suite->add(QUANTLIB_TEST_CASE(&Solver1DTest::testSecant));
223 return suite;
224}
225
226

source code of quantlib/test-suite/solvers.cpp