[go: up one dir, main page]
1/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2
3/*
4 Copyright (C) 2015 Klaus Spanderen
5
6 This file is part of QuantLib, a free-software/open-source library
7 for financial quantitative analysts and developers - http://quantlib.org/
8
9 QuantLib is free software: you can redistribute it and/or modify it
10 under the terms of the QuantLib license. You should have received a
11 copy of the license along with this program; if not, please email
12 <quantlib-dev@lists.sf.net>. The license is also available online at
13 <http://quantlib.org/license.shtml>.
14
15 This program is distributed in the hope that it will be useful, but WITHOUT
16 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
17 FOR A PARTICULAR PURPOSE. See the license for more details.
18*/
19
20#include "numericaldifferentiation.hpp"
21#include "utilities.hpp"
22
23#include <ql/math/matrix.hpp>
24#include <ql/math/factorial.hpp>
25#include <ql/methods/finitedifferences/operators/numericaldifferentiation.hpp>
26#include <cmath>
27#include <algorithm>
28
29using namespace QuantLib;
30using namespace boost::unit_test_framework;
31
32namespace {
33 bool isTheSame(Real a, Real b) {
34 constexpr double eps = 500 * QL_EPSILON;
35
36 if (std::fabs(x: b) < QL_EPSILON)
37 return std::fabs(x: a) < eps;
38 else
39 return std::fabs(x: (a - b)/b) < eps;
40 }
41
42 void checkTwoArraysAreTheSame(const Array& calculated,
43 const Array& expected) {
44 bool correct = (calculated.size() == expected.size())
45 && std::equal(first1: calculated.begin(), last1: calculated.end(),
46 first2: expected.begin(), binary_pred: isTheSame);
47
48 if (!correct) {
49 BOOST_FAIL("Failed to reproduce expected array"
50 << "\n calculated: " << calculated
51 << "\n expected: " << expected
52 << "\n difference: " << expected - calculated);
53 }
54 }
55
56 void singleValueTest(const std::string& comment,
57 Real calculated, Real expected, Real tol) {
58 if (std::fabs(x: calculated - expected) > tol)
59 BOOST_FAIL("Failed to reproduce " << comment
60 << " order derivative"
61 << "\n calculated: " << calculated
62 << "\n expected: " << expected
63 << "\n tolerance: " << tol
64 << "\n difference: "
65 << expected - calculated);
66 }
67}
68
69void NumericalDifferentiationTest::testTabulatedCentralScheme() {
70 BOOST_TEST_MESSAGE("Testing numerical differentiation "
71 "using the central scheme...");
72 const ext::function<Real(Real)> f;
73
74 const NumericalDifferentiation::Scheme central
75 = NumericalDifferentiation::Central;
76
77 // see http://en.wikipedia.org/wiki/Finite_difference_coefficient
78 checkTwoArraysAreTheSame(
79 calculated: NumericalDifferentiation(f, 1, 1.0, 3, central).weights(),
80 expected: {-0.5, 0.0, 0.5});
81
82 checkTwoArraysAreTheSame(
83 calculated: NumericalDifferentiation(f, 1, 0.5, 3, central).weights(),
84 expected: {-1.0, 0.0, 1.0});
85
86 checkTwoArraysAreTheSame(
87 calculated: NumericalDifferentiation(f, 1, 0.25, 7, central).weights(),
88 expected: {-4/60.0, 12/20.0, -12/4.0, 0.0, 12/4.0, -12/20.0, 4/60.0});
89
90 checkTwoArraysAreTheSame(
91 calculated: NumericalDifferentiation(f, 4, std::pow(x: 0.5, y: 0.25), 9, central).weights(),
92 expected: {14/240.0, -4/5.0, 338/60.0, -244/15.0, 182/8.0, -244/15.0, 338/60.0, -4/5.0, 14/240.0});
93
94 checkTwoArraysAreTheSame(
95 calculated: NumericalDifferentiation(f, 1, 0.5, 7, central).offsets(),
96 expected: {-1.5, -1.0, -0.5, 0.0, 0.5, 1.0, 1.5});
97}
98
99void NumericalDifferentiationTest::testTabulatedBackwardScheme() {
100 BOOST_TEST_MESSAGE("Testing numerical differentiation "
101 "using the backward scheme...");
102 const ext::function<Real(Real)> f;
103
104 const NumericalDifferentiation::Scheme backward
105 = NumericalDifferentiation::Backward;
106
107 // see http://en.wikipedia.org/wiki/Finite_difference_coefficient
108 checkTwoArraysAreTheSame(
109 calculated: NumericalDifferentiation(f, 1, 1.0, 2, backward).weights(),
110 expected: {1.0, -1.0});
111
112 checkTwoArraysAreTheSame(
113 calculated: NumericalDifferentiation(f, 2, 2.0, 4, backward).weights(),
114 expected: {2/4.0, -5/4.0, 4/4.0, -1.0/4.0});
115
116 checkTwoArraysAreTheSame(
117 calculated: NumericalDifferentiation(f, 4, 1.0, 6, backward).weights(),
118 expected: {3.0, -14.0, 26.0, -24.0, 11.0, -2.0});
119
120 checkTwoArraysAreTheSame(
121 calculated: NumericalDifferentiation(f, 2, 0.5, 4, backward).offsets(),
122 expected: {0.0, -0.5, -1.0, -1.5});
123}
124
125
126void NumericalDifferentiationTest::testTabulatedForwardScheme() {
127 BOOST_TEST_MESSAGE("Testing numerical differentiation "
128 "using the Forward scheme...");
129 const ext::function<Real(Real)> f;
130
131 const NumericalDifferentiation::Scheme forward
132 = NumericalDifferentiation::Forward;
133
134 // see http://en.wikipedia.org/wiki/Finite_difference_coefficient
135 checkTwoArraysAreTheSame(
136 calculated: NumericalDifferentiation(f, 1, 1.0, 2, forward).weights(),
137 expected: {-1.0, 1.0});
138
139 checkTwoArraysAreTheSame(
140 calculated: NumericalDifferentiation(f, 1, 0.5, 3, forward).weights(),
141 expected: {-6/2.0, 4.0, -2/2.0});
142
143 checkTwoArraysAreTheSame(
144 calculated: NumericalDifferentiation(f, 1, 0.5, 7, forward).weights(),
145 expected: {-98/20.0, 12.0, -30/2.0, 40/3.0, -30/4.0, 12/5.0, -2/6.0});
146
147 checkTwoArraysAreTheSame(
148 calculated: NumericalDifferentiation(f, 2, 0.5, 4, forward).offsets(),
149 expected: {0.0, 0.5, 1.0, 1.5});
150}
151
152
153void NumericalDifferentiationTest::testIrregularSchemeFirstOrder() {
154 BOOST_TEST_MESSAGE("Testing numerical differentiation "
155 "of first order using an irregular scheme...");
156 const ext::function<Real(Real)> f;
157
158 const Real h1 = 5e-7;
159 const Real h2 = 3e-6;
160
161 const Real alpha = -h2/(h1*(h1+h2));
162 const Real gamma = h1/(h2*(h1+h2));
163 const Real beta = -alpha - gamma;
164
165 Array offsets = { -h1, 0.0, h2 };
166
167 checkTwoArraysAreTheSame(
168 calculated: NumericalDifferentiation(f, 1, offsets).weights(),
169 expected: { alpha, beta, gamma });
170}
171
172void NumericalDifferentiationTest::testIrregularSchemeSecondOrder() {
173 BOOST_TEST_MESSAGE("Testing numerical differentiation "
174 "of second order using an irregular scheme...");
175 const ext::function<Real(Real)> f;
176
177 const Real h1 = 2e-7;
178 const Real h2 = 8e-8;
179
180 const Real alpha = 2/(h1*(h1+h2));
181 const Real gamma = 2/(h2*(h1+h2));
182 const Real beta = -alpha - gamma;
183
184 Array offsets = { -h1, 0.0, h2 };
185
186 checkTwoArraysAreTheSame(
187 calculated: NumericalDifferentiation(f, 2, offsets).weights(),
188 expected: {alpha, beta, gamma});
189}
190
191
192void NumericalDifferentiationTest::testDerivativesOfSineFunction() {
193 BOOST_TEST_MESSAGE("Testing numerical differentiation"
194 " of sin function...");
195
196 const ext::function<Real(Real)> f = [](Real x) -> Real { return std::sin(x: x); };
197
198 const ext::function<Real(Real)> df_central
199 = NumericalDifferentiation(f, 1, std::sqrt(QL_EPSILON), 3,
200 NumericalDifferentiation::Central);
201
202 const ext::function<Real(Real)> df_backward
203 = NumericalDifferentiation(f, 1, std::sqrt(QL_EPSILON), 3,
204 NumericalDifferentiation::Backward);
205
206 const ext::function<Real(Real)> df_forward
207 = NumericalDifferentiation(f, 1, std::sqrt(QL_EPSILON), 3,
208 NumericalDifferentiation::Forward);
209
210 for (Real x=0.0; x < 5.0; x+=0.1) {
211 const Real calculatedCentral = df_central(x);
212 const Real calculatedBackward = df_backward(x);
213 const Real calculatedForward = df_forward(x);
214 const Real expected = std::cos(x: x);
215
216 singleValueTest(comment: "central first", calculated: calculatedCentral, expected, tol: 1e-8);
217 singleValueTest(comment: "backward first", calculated: calculatedBackward, expected, tol: 1e-6);
218 singleValueTest(comment: "forward first", calculated: calculatedForward, expected, tol: 1e-6);
219 }
220
221 const ext::function<Real(Real)> df4_central
222 = NumericalDifferentiation(f, 4, 1e-2, 7,
223 NumericalDifferentiation::Central);
224 const ext::function<Real(Real)> df4_backward
225 = NumericalDifferentiation(f, 4, 1e-2, 7,
226 NumericalDifferentiation::Backward);
227 const ext::function<Real(Real)> df4_forward
228 = NumericalDifferentiation(f, 4, 1e-2, 7,
229 NumericalDifferentiation::Forward);
230
231 for (Real x=0.0; x < 5.0; x+=0.1) {
232 const Real calculatedCentral = df4_central(x);
233 const Real calculatedBackward = df4_backward(x);
234 const Real calculatedForward = df4_forward(x);
235 const Real expected = std::sin(x: x);
236
237 singleValueTest(comment: "central 4th", calculated: calculatedCentral, expected, tol: 1e-4);
238 singleValueTest(comment: "backward 4th", calculated: calculatedBackward, expected, tol: 1e-4);
239 singleValueTest(comment: "forward 4th", calculated: calculatedForward, expected, tol: 1e-4);
240 }
241
242 const Array offsets = {-0.01, -0.02, 0.03, 0.014, 0.041};
243 NumericalDifferentiation df3_irregular(f, 3, offsets);
244
245 checkTwoArraysAreTheSame(calculated: df3_irregular.offsets(), expected: offsets);
246
247 for (Real x=0.0; x < 5.0; x+=0.1) {
248 const Real calculatedIrregular = df3_irregular(x);
249 const Real expected = -std::cos(x: x);
250
251 singleValueTest(comment: "irregular 3th", calculated: calculatedIrregular, expected, tol: 5e-5);
252 }
253}
254
255namespace {
256 Array vandermondeCoefficients(
257 Size order, Real x, const Array& gridPoints) {
258
259 const Array q = gridPoints - x;
260 const Size n = gridPoints.size();
261
262 Matrix m(n, n, 1.0);
263 for (Size i=1; i < n; ++i) {
264 const Real fact = Factorial::get(n: i);
265 for (Size j=0; j < n; ++j)
266 m[i][j] = std::pow(x: q[j], y: Integer(i)) / fact;
267 }
268
269 Array b(n, 0.0);
270 b[order] = 1.0;
271 return inverse(m)*b;
272 }
273}
274
275void NumericalDifferentiationTest::testCoefficientBasedOnVandermonde() {
276 BOOST_TEST_MESSAGE("Testing coefficients from numerical differentiation"
277 " by comparison with results from"
278 " Vandermonde matrix inversion...");
279 const ext::function<Real(Real)> f;
280
281 for (Natural order=0; order < 5; ++order) {
282 for (Natural nGridPoints = order + 1;
283 nGridPoints < order + 3; ++nGridPoints) {
284
285 Array gridPoints(nGridPoints);
286 for (Natural i=0; i < nGridPoints; ++i) {
287 const Real p = Real(i);
288 gridPoints[i] = std::sin(x: p) + std::cos(x: p); // strange points
289 }
290
291 const Real x = 0.3902842; // strange points
292 const Array weightsVandermonde
293 = vandermondeCoefficients(order, x, gridPoints);
294 const NumericalDifferentiation nd(f, order, gridPoints-x);
295
296 checkTwoArraysAreTheSame(calculated: gridPoints, expected: nd.offsets() + x);
297 checkTwoArraysAreTheSame(calculated: weightsVandermonde, expected: nd.weights());
298 }
299 }
300}
301
302
303
304test_suite* NumericalDifferentiationTest::suite() {
305 auto* suite = BOOST_TEST_SUITE("NumericalDifferentiation tests");
306
307 suite->add(QUANTLIB_TEST_CASE(
308 &NumericalDifferentiationTest::testTabulatedCentralScheme));
309 suite->add(QUANTLIB_TEST_CASE(
310 &NumericalDifferentiationTest::testTabulatedBackwardScheme));
311 suite->add(QUANTLIB_TEST_CASE(
312 &NumericalDifferentiationTest::testTabulatedForwardScheme));
313 suite->add(QUANTLIB_TEST_CASE(
314 &NumericalDifferentiationTest::testIrregularSchemeFirstOrder));
315 suite->add(QUANTLIB_TEST_CASE(
316 &NumericalDifferentiationTest::testIrregularSchemeSecondOrder));
317 suite->add(QUANTLIB_TEST_CASE(
318 &NumericalDifferentiationTest::testDerivativesOfSineFunction));
319 suite->add(QUANTLIB_TEST_CASE(
320 &NumericalDifferentiationTest::testCoefficientBasedOnVandermonde));
321
322 return suite;
323}
324
325
326

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