QuantLib: a free/open-source library for quantitative finance
fully annotated source code - version 1.34
Loading...
Searching...
No Matches
concentrating1dmesher.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) 2009 Ralph Schreyer
5 Copyright (C) 2014 Johannes Göttker-Schnetmann
6 Copyright (C) 2014 Klaus Spanderen
7 Copyright (C) 2015 Peter Caspers
8
9 This file is part of QuantLib, a free-software/open-source library
10 for financial quantitative analysts and developers - http://quantlib.org/
11
12 QuantLib is free software: you can redistribute it and/or modify it
13 under the terms of the QuantLib license. You should have received a
14 copy of the license along with this program; if not, please email
15 <quantlib-dev@lists.sf.net>. The license is also available online at
16 <http://quantlib.org/license.shtml>.
17
18 This program is distributed in the hope that it will be useful, but WITHOUT
19 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
20 FOR A PARTICULAR PURPOSE. See the license for more details.
21*/
22
23/*! \file concentrating1dmesher.cpp
24 \brief One-dimensional grid mesher concentrating around critical points
25*/
26
27#include <ql/errors.hpp>
28#include <ql/timegrid.hpp>
29#include <ql/utilities/null.hpp>
30#include <ql/math/array.hpp>
37#include <cmath>
38
39namespace QuantLib {
40
42 Real start, Real end, Size size, const std::pair<Real, Real>& cPoints,
43 const bool requireCPoint)
44 : Fdm1dMesher(size) {
45
46 QL_REQUIRE(end > start, "end must be larger than start");
47
48 const Real cPoint = cPoints.first;
49 const Real density = cPoints.second == Null<Real>() ?
50 Null<Real>() : cPoints.second*(end - start);
51
52 QL_REQUIRE(cPoint == Null<Real>() || (cPoint >= start && cPoint <= end),
53 "cPoint must be between start and end");
54 QL_REQUIRE(density == Null<Real>() || density > 0.0,
55 "density > 0 required");
56 QL_REQUIRE(cPoint == Null<Real>() || density != Null<Real>(),
57 "density must be given if cPoint is given");
58 QL_REQUIRE(!requireCPoint || cPoint != Null<Real>(),
59 "cPoint is required in grid but not given");
60
61 const Real dx = 1.0 / (size - 1);
62
63 if (cPoint != Null<Real>()) {
64 std::vector<Real> u, z;
65 ext::shared_ptr<Interpolation> transform;
66 const Real c1 = std::asinh((start - cPoint) / density);
67 const Real c2 = std::asinh((end - cPoint) / density);
68 if (requireCPoint) {
69 u.push_back(0.0);
70 z.push_back(0.0);
71 if (!close(cPoint, start) && !close(cPoint, end)) {
72 const Real z0 = -c1 / (c2 - c1);
73 const Real u0 =
74 std::max(std::min(std::lround(z0 * (size - 1)),
75 static_cast<long>(size) - 2),
76 1L) /
77 ((Real)(size - 1));
78 u.push_back(u0);
79 z.push_back(z0);
80 }
81 u.push_back(1.0);
82 z.push_back(1.0);
83 transform = ext::shared_ptr<Interpolation>(
84 new LinearInterpolation(u.begin(), u.end(), z.begin()));
85 }
86
87 for (Size i = 1; i < size - 1; ++i) {
88 const Real li = requireCPoint ? (*transform)(i*dx) : i*dx;
89 locations_[i] = cPoint
90 + density*std::sinh(c1*(1.0 - li) + c2*li);
91 }
92 }
93 else {
94 for (Size i = 1; i < size - 1; ++i) {
95 locations_[i] = start + i*dx*(end - start);
96 }
97 }
98
99 locations_.front() = start;
100 locations_.back() = end;
101
102 for (Size i = 0; i < size - 1; ++i) {
103 dplus_[i] = dminus_[i + 1] = locations_[i + 1] - locations_[i];
104 }
105 dplus_.back() = dminus_.front() = Null<Real>();
106 }
107
108 namespace {
109 class OdeIntegrationFct {
110 public:
111 OdeIntegrationFct(const std::vector<Real>& points,
112 const std::vector<Real>& betas,
113 Real tol)
114 : rk_(tol), points_(points), betas_(betas) {}
115
116 Real solve(Real a, Real y0, Real x0, Real x1) {
117 AdaptiveRungeKutta<>::OdeFct1d odeFct([&](Real x, Real y){ return jac(a, x, y); });
118 return rk_(odeFct, y0, x0, x1);
119 }
120
121 private:
122 Real jac(Real a, Real, Real y) const {
123 Real s=0.0;
124 for (Size i=0; i < points_.size(); ++i) {
125 s+=1.0/(betas_[i] + squared(y - points_[i]));
126 }
127 return a/std::sqrt(s);
128 }
129
130 AdaptiveRungeKutta<> rk_;
131 const std::vector<Real> &points_, &betas_;
132 };
133
134 bool equal_on_first(const std::pair<Real, Real>& p1,
135 const std::pair<Real, Real>& p2) {
136 return close_enough(p1.first, p2.first, 1000);
137 }
138 }
139
140
142 Real start, Real end, Size size,
143 const std::vector<ext::tuple<Real, Real, bool> >& cPoints,
144 Real tol)
145 : Fdm1dMesher(size) {
146
147 QL_REQUIRE(end > start, "end must be larger than start");
148
149 std::vector<Real> points, betas;
150 for (const auto& cPoint : cPoints) {
151 points.push_back(ext::get<0>(cPoint));
152 betas.push_back(squared(ext::get<1>(cPoint) * (end - start)));
153 }
154
155 // get scaling factor a so that y(1) = end
156 Real aInit = 0.0;
157 for (Size i=0; i < points.size(); ++i) {
158 const Real c1 = std::asinh((start-points[i])/betas[i]);
159 const Real c2 = std::asinh((end-points[i])/betas[i]);
160 aInit+=(c2-c1)/points.size();
161 }
162
163 OdeIntegrationFct fct(points, betas, tol);
164 const Real a = Brent().solve(
165 [&](Real x) -> Real { return fct.solve(x, start, 0.0, 1.0) - end; },
166 tol, aInit, 0.1*aInit);
167
168 // solve ODE for all grid points
169 Array x(size), y(size);
170 x[0] = 0.0; y[0] = start;
171 const Real dx = 1.0/(size-1);
172 for (Size i=1; i < size; ++i) {
173 x[i] = i*dx;
174 y[i] = fct.solve(a, y[i-1], x[i-1], x[i]);
175 }
176
177 // eliminate numerical noise and ensure y(1) = end
178 const Real dy = y.back() - end;
179 for (Size i=1; i < size; ++i) {
180 y[i] -= i*dx*dy;
181 }
182
183 LinearInterpolation odeSolution(x.begin(), x.end(), y.begin());
184
185 // ensure required points are part of the grid
186 std::vector<std::pair<Real, Real> > w(1, std::make_pair(0.0, 0.0));
187
188 for (Size i=0; i < points.size(); ++i) {
189 if (ext::get<2>(cPoints[i]) && points[i] > start && points[i] < end) {
190
191 const Size j = std::distance(y.begin(),
192 std::lower_bound(y.begin(), y.end(), points[i]));
193
194 const Real e = Brent().solve(
195 [&](Real x) -> Real { return odeSolution(x, true) - points[i]; },
196 QL_EPSILON, x[j], 0.5/size);
197
198 w.emplace_back(std::min(x[size - 2], x[j]), e);
199 }
200 }
201 w.emplace_back(1.0, 1.0);
202 std::sort(w.begin(), w.end());
203 w.erase(std::unique(w.begin(), w.end(), equal_on_first), w.end());
204
205 std::vector<Real> u(w.size()), z(w.size());
206 for (Size i=0; i < w.size(); ++i) {
207 u[i] = w[i].first;
208 z[i] = w[i].second;
209 }
210 LinearInterpolation transform(u.begin(), u.end(), z.begin());
211
212 for (Size i=0; i < size; ++i) {
213 locations_[i] = odeSolution(transform(i*dx));
214 }
215
216 for (Size i=0; i < size-1; ++i) {
217 dplus_[i] = dminus_[i+1] = locations_[i+1] - locations_[i];
218 }
219 dplus_.back() = dminus_.front() = Null<Real>();
220 }
221}
Runge-Kutta ODE integration.
1-D array used in linear algebra.
Brent 1-D solver.
ext::function< T(const Real, const T)> OdeFct1d
1-D array used in linear algebra.
Definition: array.hpp:52
Brent 1-D solver
Definition: brent.hpp:37
Concentrating1dMesher(Real start, Real end, Size size, const std::pair< Real, Real > &cPoints=(std::pair< Real, Real >(Null< Real >(), Null< Real >())), bool requireCPoint=false)
std::vector< Real > locations_
Definition: fdm1dmesher.hpp:47
std::vector< Real > dplus_
Definition: fdm1dmesher.hpp:48
std::vector< Real > dminus_
Definition: fdm1dmesher.hpp:48
Linear interpolation between discrete points
template class providing a null value for a given type.
Definition: null.hpp:76
Real solve(const F &f, Real accuracy, Real guess, Real step) const
Definition: solver1d.hpp:84
floating-point comparisons
const std::vector< Real > & betas_
const std::vector< Real > & points_
AdaptiveRungeKutta rk_
One-dimensional grid mesher concentrating around critical points.
Classes and functions for error handling.
#define QL_REQUIRE(condition, message)
throw an error if the given pre-condition is not verified
Definition: errors.hpp:117
#define QL_EPSILON
Definition: qldefines.hpp:178
QL_REAL Real
real number
Definition: types.hpp:50
std::size_t Size
size of a container
Definition: types.hpp:58
linear interpolation between discrete points
functionals and combinators not included in the STL
Definition: any.hpp:35
T squared(T x)
Definition: functional.hpp:37
bool close(const Quantity &m1, const Quantity &m2, Size n)
Definition: quantity.cpp:163
bool close_enough(const Quantity &m1, const Quantity &m2, Size n)
Definition: quantity.cpp:182
null values
discrete time grid