QuantLib: a free/open-source library for quantitative finance
fully annotated source code - version 1.34
Loading...
Searching...
No Matches
numericaldifferentiation.cpp
Go to the documentation of this file.
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/*! \file numericaldifferentiation.cpp */
21
23
24#pragma push_macro("BOOST_DISABLE_ASSERTS")
25#ifndef BOOST_DISABLE_ASSERTS
26#define BOOST_DISABLE_ASSERTS
27#endif
28#include <boost/multi_array.hpp>
29#pragma pop_macro("BOOST_DISABLE_ASSERTS")
30
31#include <utility>
32
33namespace QuantLib {
34
35 namespace {
36 Array calcOffsets(
38 QL_REQUIRE(n > 1, "number of steps must be greater than one");
39
40 Array retVal(n);
41 switch (scheme) {
43 QL_REQUIRE(n > 2 && (n % 2),
44 "number of steps must be an odd number greater than two");
45 for (Integer i=0; i < Integer(n); ++i)
46 retVal[i] = (i-Integer(n/2))*h;
47 break;
49 for (Size i=0; i < n; ++i)
50 retVal[i]=-(i*h);
51 break;
53 for (Size i=0; i < n; ++i)
54 retVal[i]=i*h;
55 break;
56 default:
57 QL_FAIL("unknown numerical differentiation scheme");
58 }
59
60 return retVal;
61 }
62
63 // This is a C++ implementation of the algorithm/pseudo code in
64 // B. Fornberg, 1998. Calculation of Weights
65 // in Finite Difference Formulas
66 // https://amath.colorado.edu/faculty/fornberg/Docs/sirev_cl.pdf
67 Array calcWeights(const Array& x, Size M) {
68 const Size N = x.size();
69 QL_REQUIRE(N > M, "number of points must be greater "
70 "than the order of the derivative");
71
72 boost::multi_array<Real, 3> d(boost::extents[M+1][N][N]);
73 d[0][0][0] = 1.0;
74 Real c1 = 1.0;
75
76 for (Size n=1; n < N; ++n) {
77 Real c2 = 1.0;
78 for (Size nu=0; nu < n; ++nu) {
79 const Real c3 = x[n] - x[nu];
80 c2*=c3;
81
82 for (Size m=0; m <= std::min(n, M); ++m) {
83 d[m][n][nu] = (x[n]*d[m][n-1][nu]
84 - ((m > 0)? Real(m*d[m-1][n-1][nu]) : 0.0))/c3;
85 }
86 }
87
88 for (Size m=0; m <= M; ++m) {
89 d[m][n][n] = c1/c2*( ((m > 0)? Real(m*d[m-1][n-1][n-1]) : 0.0) -
90 x[n-1]*d[m][n-1][n-1] );
91 }
92 c1=c2;
93 }
94
95
96 Array retVal(N);
97 for (Size i=0; i < N; ++i) {
98 retVal[i] = d[M][N-1][i];
99 }
100 return retVal;
101 }
102 }
103
105 Size orderOfDerivative,
106 Array x_offsets)
107 : offsets_(std::move(x_offsets)), w_(calcWeights(offsets_, orderOfDerivative)),
108 f_(std::move(f)) {}
109
110
112 Size orderOfDerivative,
113 Real stepSize,
114 Size steps,
115 Scheme scheme)
116 : offsets_(calcOffsets(stepSize, steps, scheme)), w_(calcWeights(offsets_, orderOfDerivative)),
117 f_(std::move(f)) {}
118}
1-D array used in linear algebra.
Definition: array.hpp:52
NumericalDifferentiation(ext::function< Real(Real)> f, Size orderOfDerivative, Array x_offsets)
ext::function< Real(Real)> f_
#define QL_REQUIRE(condition, message)
throw an error if the given pre-condition is not verified
Definition: errors.hpp:117
#define QL_FAIL(message)
throw an error (possibly with file and line information)
Definition: errors.hpp:92
Date d
QL_REAL Real
real number
Definition: types.hpp:50
QL_INTEGER Integer
integer number
Definition: types.hpp:35
std::size_t Size
size of a container
Definition: types.hpp:58
Definition: any.hpp:35
STL namespace.
numerical differentiation of arbitrary order and on irregular grids
Real nu
Definition: sabr.cpp:200