24#ifndef quantlib_tr_bdf2_scheme_hpp
25#define quantlib_tr_bdf2_scheme_hpp
38 template <
class Trapezo
idalScheme>
52 ext::shared_ptr<FdmLinearOpComposite> map,
53 const ext::shared_ptr<TrapezoidalScheme>& trapezoidalScheme,
70 const ext::shared_ptr<FdmLinearOpComposite>
map_;
77 template <
class Trapezo
idalScheme>
80 ext::shared_ptr<FdmLinearOpComposite> map,
81 const ext::shared_ptr<TrapezoidalScheme>& trapezoidalScheme,
86 alpha_(
alpha), map_(
std::move(map)), trapezoidalScheme_(trapezoidalScheme), bcSet_(bcSet),
87 relTol_(relTol), solverType_(solverType) {}
89 template <
class Trapezo
idalScheme>
92 beta_= (1.0-alpha_)/(2.0-alpha_)*dt_;
95 template <
class Trapezo
idalScheme>
100 template <
class Trapezo
idalScheme>
102 return r - beta_*map_->apply(
r);
105 template <
class Trapezo
idalScheme>
107 QL_REQUIRE(
t-dt_ > -1e-8,
"a step towards negative time given");
109 const Time intermediateTimeStep = dt_*alpha_;
112 trapezoidalScheme_->setStep(intermediateTimeStep);
113 trapezoidalScheme_->step(fStar,
t);
115 bcSet_.setTime(std::max(0.0,
t-dt_));
116 bcSet_.applyBeforeSolving(*map_, fn);
119 (1/alpha_*fStar -
squared(1-alpha_)/alpha_*fn)/(2-alpha_);
121 if (map_->size() == 1) {
122 fn = map_->solve_splitting(0,
f, -beta_);
125 auto preconditioner = [&](
const Array& _a){
return map_->preconditioner(_a, -beta_); };
126 auto applyF = [&](
const Array& _a){
return apply(_a); };
131 relTol_, preconditioner).
solve(
f,
f);
135 }
else if (solverType_ ==
GMRES) {
141 (*iterations_) += result.
errors.size();
145 QL_FAIL(
"unknown/illegal solver type");
148 bcSet_.applyAfterSolving(fn);
Biconjugate gradient stabilized method.
1-D array used in linear algebra.
BiCGStabResult solve(const Array &b, const Array &x0=Array()) const
GMRESResult solve(const Array &b, const Array &x0=Array()) const
template class providing a null value for a given type.
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_
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)
#define QL_REQUIRE(condition, message)
throw an error if the given pre-condition is not verified
#define QL_FAIL(message)
throw an error (possibly with file and line information)
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
std::size_t Size
size of a container
functionals and combinators not included in the STL
Differential operator traits.
ext::shared_ptr< YieldTermStructure > r