24#ifndef quantlib_tr_bdf2_scheme_hpp
25#define quantlib_tr_bdf2_scheme_hpp
27#include <ql/functional.hpp>
28#include <ql/math/functional.hpp>
29#include <ql/math/matrixutilities/bicgstab.hpp>
30#include <ql/math/matrixutilities/gmres.hpp>
31#include <ql/methods/finitedifferences/operators/fdmlinearopcomposite.hpp>
32#include <ql/methods/finitedifferences/operatortraits.hpp>
33#include <ql/methods/finitedifferences/schemes/boundaryconditionschemehelper.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);
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)
Real Time
continuous quantity with 1-year units
std::size_t Size
size of a container