| 1 | /* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */ |
|---|
| 2 | |
|---|
| 3 | /* |
|---|
| 4 | Copyright (C) 2014 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 <ql/math/integrals/discreteintegrals.hpp> |
|---|
| 21 | #include <boost/accumulators/accumulators.hpp> |
|---|
| 22 | #include <boost/accumulators/statistics/sum.hpp> |
|---|
| 23 | |
|---|
| 24 | using namespace boost::accumulators; |
|---|
| 25 | |
|---|
| 26 | namespace QuantLib { |
|---|
| 27 | |
|---|
| 28 | Real DiscreteTrapezoidIntegral::operator()( |
|---|
| 29 | const Array& x, const Array& f) const { |
|---|
| 30 | |
|---|
| 31 | const Size n = f.size(); |
|---|
| 32 | QL_REQUIRE(n == x.size(), "inconsistent size" ); |
|---|
| 33 | |
|---|
| 34 | accumulator_set<Real, features<tag::sum> > acc; |
|---|
| 35 | |
|---|
| 36 | for (Size i=0; i < n-1; ++i) { |
|---|
| 37 | acc((x[i+1]-x[i])*(f[i]+f[i+1])); |
|---|
| 38 | } |
|---|
| 39 | |
|---|
| 40 | return 0.5*sum(acc); |
|---|
| 41 | } |
|---|
| 42 | |
|---|
| 43 | Real DiscreteSimpsonIntegral::operator()( |
|---|
| 44 | const Array& x, const Array& f) const { |
|---|
| 45 | |
|---|
| 46 | const Size n = f.size(); |
|---|
| 47 | QL_REQUIRE(n == x.size(), "inconsistent size" ); |
|---|
| 48 | |
|---|
| 49 | accumulator_set<Real, features<tag::sum> > acc; |
|---|
| 50 | |
|---|
| 51 | for (Size j=0; j < n-2; j+=2) { |
|---|
| 52 | const Real dxj = x[j+1]-x[j]; |
|---|
| 53 | const Real dxjp1 = x[j+2]-x[j+1]; |
|---|
| 54 | |
|---|
| 55 | const Real alpha = -dxjp1*(2*x[j]-3*x[j+1]+x[j+2]); |
|---|
| 56 | const Real dd = x[j+2]-x[j]; |
|---|
| 57 | const Real k = dd/(6*dxjp1*dxj); |
|---|
| 58 | const Real beta = dd*dd; |
|---|
| 59 | const Real gamma = dxj*(x[j]-3*x[j+1]+2*x[j+2]); |
|---|
| 60 | |
|---|
| 61 | acc(k*alpha*f[j]+k*beta*f[j+1]+k*gamma*f[j+2]); |
|---|
| 62 | } |
|---|
| 63 | if ((n & 1) == 0U) { |
|---|
| 64 | acc(0.5*(x[n-1]-x[n-2])*(f[n-1]+f[n-2])); |
|---|
| 65 | } |
|---|
| 66 | |
|---|
| 67 | return sum(acc); |
|---|
| 68 | } |
|---|
| 69 | |
|---|
| 70 | |
|---|
| 71 | Real DiscreteTrapezoidIntegrator::integrate( |
|---|
| 72 | const ext::function<Real (Real)>& f, Real a, Real b) const { |
|---|
| 73 | const Array x(maxEvaluations(), a, (b-a)/(maxEvaluations()-1)); |
|---|
| 74 | Array fv(x.size()); |
|---|
| 75 | std::transform(first: x.begin(), last: x.end(), result: fv.begin(), unary_op: f); |
|---|
| 76 | |
|---|
| 77 | increaseNumberOfEvaluations(increase: maxEvaluations()); |
|---|
| 78 | return DiscreteTrapezoidIntegral()(x, fv); |
|---|
| 79 | } |
|---|
| 80 | |
|---|
| 81 | Real DiscreteSimpsonIntegrator::integrate( |
|---|
| 82 | const ext::function<Real (Real)>& f, Real a, Real b) const { |
|---|
| 83 | const Array x(maxEvaluations(), a, (b-a)/(maxEvaluations()-1)); |
|---|
| 84 | Array fv(x.size()); |
|---|
| 85 | std::transform(first: x.begin(), last: x.end(), result: fv.begin(), unary_op: f); |
|---|
| 86 | |
|---|
| 87 | increaseNumberOfEvaluations(increase: maxEvaluations()); |
|---|
| 88 | return DiscreteSimpsonIntegral()(x, fv); |
|---|
| 89 | } |
|---|
| 90 | } |
|---|
| 91 | |
|---|
