QuantLib: a free/open-source library for quantitative finance
fully annotated source code - version 1.34
Loading...
Searching...
No Matches
trbdf2scheme.hpp
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) 2018 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 trbdf2scheme.hpp
21 \brief trapezoidal BDF2 scheme
22*/
23
24#ifndef quantlib_tr_bdf2_scheme_hpp
25#define quantlib_tr_bdf2_scheme_hpp
26
27#include <ql/functional.hpp>
34#include <utility>
35
36namespace QuantLib {
37
38 template <class TrapezoidalScheme>
40 public:
42
43 // typedefs
49
50 // constructors
52 ext::shared_ptr<FdmLinearOpComposite> map,
53 const ext::shared_ptr<TrapezoidalScheme>& trapezoidalScheme,
54 const bc_set& bcSet = bc_set(),
55 Real relTol = 1e-8,
56 SolverType solverType = BiCGstab);
57
58 void step(array_type& a, Time t);
59 void setStep(Time dt);
60
62 protected:
63 Array apply(const Array& r) const;
64
67 ext::shared_ptr<Size> iterations_;
68
69 const Real alpha_;
70 const ext::shared_ptr<FdmLinearOpComposite> map_;
71 const ext::shared_ptr<TrapezoidalScheme>& trapezoidalScheme_;
75 };
76
77 template <class TrapezoidalScheme>
79 Real alpha,
80 ext::shared_ptr<FdmLinearOpComposite> map,
81 const ext::shared_ptr<TrapezoidalScheme>& trapezoidalScheme,
82 const bc_set& bcSet,
83 Real relTol,
84 SolverType solverType)
85 : dt_(Null<Real>()), beta_(Null<Real>()), iterations_(ext::make_shared<Size>(0U)),
86 alpha_(alpha), map_(std::move(map)), trapezoidalScheme_(trapezoidalScheme), bcSet_(bcSet),
87 relTol_(relTol), solverType_(solverType) {}
88
89 template <class TrapezoidalScheme>
91 dt_=dt;
92 beta_= (1.0-alpha_)/(2.0-alpha_)*dt_;
93 }
94
95 template <class TrapezoidalScheme>
97 return *iterations_;
98 }
99
100 template <class TrapezoidalScheme>
102 return r - beta_*map_->apply(r);
103 }
104
105 template <class TrapezoidalScheme>
107 QL_REQUIRE(t-dt_ > -1e-8, "a step towards negative time given");
108
109 const Time intermediateTimeStep = dt_*alpha_;
110
111 array_type fStar = fn;
112 trapezoidalScheme_->setStep(intermediateTimeStep);
113 trapezoidalScheme_->step(fStar, t);
114
115 bcSet_.setTime(std::max(0.0, t-dt_));
116 bcSet_.applyBeforeSolving(*map_, fn);
117
118 const array_type f =
119 (1/alpha_*fStar - squared(1-alpha_)/alpha_*fn)/(2-alpha_);
120
121 if (map_->size() == 1) {
122 fn = map_->solve_splitting(0, f, -beta_);
123 }
124 else {
125 auto preconditioner = [&](const Array& _a){ return map_->preconditioner(_a, -beta_); };
126 auto applyF = [&](const Array& _a){ return apply(_a); };
127
128 if (solverType_ == BiCGstab) {
129 const BiCGStabResult result =
130 QuantLib::BiCGstab(applyF, std::max(Size(10), fn.size()),
131 relTol_, preconditioner).solve(f, f);
132
133 (*iterations_) += result.iterations;
134 fn = result.x;
135 } else if (solverType_ == GMRES) {
136 const GMRESResult result =
137 QuantLib::GMRES(applyF, std::max(Size(10), fn.size() / 10U), relTol_,
138 preconditioner)
139 .solve(f, f);
140
141 (*iterations_) += result.errors.size();
142 fn = result.x;
143 }
144 else
145 QL_FAIL("unknown/illegal solver type");
146 }
147
148 bcSet_.applyAfterSolving(fn);
149 }
150}
151
152#endif
Biconjugate gradient stabilized method.
1-D array used in linear algebra.
Definition: array.hpp:52
BiCGStabResult solve(const Array &b, const Array &x0=Array()) const
Definition: bicgstab.cpp:37
GMRESResult solve(const Array &b, const Array &x0=Array()) const
Definition: gmres.cpp:39
template class providing a null value for a given type.
Definition: null.hpp:76
std::vector< ext::shared_ptr< bc_type > > bc_set
Operator::array_type array_type
condition to be applied at every time step
traits::operator_type operator_type
Array apply(const Array &r) const
ext::shared_ptr< Size > iterations_
const ext::shared_ptr< TrapezoidalScheme > & trapezoidalScheme_
const SolverType solverType_
const BoundaryConditionSchemeHelper bcSet_
traits::bc_set bc_set
const ext::shared_ptr< FdmLinearOpComposite > map_
traits::condition_type condition_type
traits::array_type array_type
OperatorTraits< FdmLinearOp > traits
Size numberOfIterations() const
TrBDF2Scheme(Real alpha, ext::shared_ptr< FdmLinearOpComposite > map, const ext::shared_ptr< TrapezoidalScheme > &trapezoidalScheme, const bc_set &bcSet=bc_set(), Real relTol=1e-8, SolverType solverType=BiCGstab)
void step(array_type &a, Time t)
const DefaultType & t
#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
composite pattern for linear operators
Maps function, bind and cref to either the boost or std implementation.
generalized minimal residual method
Real Time
continuous quantity with 1-year units
Definition: types.hpp:62
QL_REAL Real
real number
Definition: types.hpp:50
std::size_t Size
size of a container
Definition: types.hpp:58
functionals and combinators not included in the STL
Definition: any.hpp:35
T squared(T x)
Definition: functional.hpp:37
STL namespace.
Differential operator traits.
ext::shared_ptr< YieldTermStructure > r
Real alpha
Definition: sabr.cpp:200
std::list< Real > errors
Definition: gmres.hpp:50