26#include <ql/methods/finitedifferences/meshers/fdmmesher.hpp>
27#include <ql/methods/finitedifferences/operators/fdmlinearoplayout.hpp>
28#include <ql/methods/finitedifferences/operators/firstderivativeop.hpp>
29#include <ql/methods/finitedifferences/operators/secondderivativeop.hpp>
30#include <ql/methods/finitedifferences/operators/fdmsquarerootfwdop.hpp>
31#include <ql/methods/finitedifferences/operators/modtriplebandlinearop.hpp>
36 const ext::shared_ptr<FdmMesher>& mesher,
39 : direction_(direction),
43 transform_(transform),
44 mapX_(transform == Plain ?
46 .mult(kappa*(mesher->locations(direction_)-theta) + sigma*sigma)
48 .mult(0.5*sigma*sigma*mesher->locations(direction_)))
49 .add(
Array(mesher->layout()->size(), kappa)))
53 .mult(0.5*sigma*sigma*mesher->locations(direction_))
55 .mult(kappa*(mesher->locations(direction_)+theta)))
56 .add(
Array(mesher->layout()->size(),
57 2*kappa*kappa*theta/(sigma*sigma))))
60 .mult(
Exp(-mesher->locations(direction))
61 *( -0.5*sigma*sigma - kappa*theta) + kappa)
63 .mult(0.5*sigma*sigma*
Exp(-mesher->locations(direction))))
64 .add(kappa*theta*
Exp(-mesher->locations(direction))))
66 v_ (mesher->layout()->dim()[direction_]) {
68 for (
const auto& iter : *mesher->layout()) {
79 const ext::shared_ptr<FdmMesher>& mesher) {
81 Real alpha, beta, gamma;
89 for (
const auto& iter : *mesher->layout()) {
91 const Size idx = iter.index();
92 mapX_->diag(idx) = beta + f*b;
93 mapX_->upper(idx) = gamma + f*c;
99 const ext::shared_ptr<FdmMesher>& mesher) {
101 Real alpha, beta, gamma;
109 for (
const auto& iter : *mesher->layout()) {
111 const Size idx = iter.index();
112 mapX_->diag(idx) = beta + f*b;
113 mapX_->lower(idx) = alpha + f*c;
119 if (transform ==
Plain) {
122 else if (transform ==
Power) {
125 else if (transform ==
Log) {
129 QL_FAIL(
"unknown transform");
133 if (transform ==
Plain) {
136 else if (transform ==
Power) {
139 else if (transform ==
Log) {
143 QL_FAIL(
"unknown transform");
153 return alpha/nu*
v(n-1);
163 return gamma/nu*
v(n+1);
173 return alpha/nu*
v(n-1);
183 return gamma/nu*
v(n+1);
194 return alpha/nu*exp(-
v(n-1));
205 return gamma/nu*exp(-
v(n+1));
209 if (i > 0 && i <=
v_.
size()) {
214 return 2*
v_[0] -
v_[1];
217 return std::max(0.5*
v_[0],
v_[0] - 0.01 * (
v_[1] -
v_[0]));
220 else if (i ==
v_.
size()+1) {
224 QL_FAIL(
"unknown index");
229 return v(i+1) -
v(i);
235 return h(i-1)*(
h(i-1)+
h(i));
241 return h(i)*(
h(i-1)+
h(i));
291 return mapX_->apply(p);
295 return Array(r.size(), 0.0);
301 return mapX_->apply(r);
304 return Array(r.size(), 0.0);
311 return mapX_->solve_splitting(r, dt, 1.0);
324 return std::vector<SparseMatrix>(1,
mapX_->toMatrix());
1-D array used in linear algebra.
const_iterator end() const
Size size() const
dimension of the array
Size size() const override
Array apply_direction(Size direction, const Array &r) const override
Array preconditioner(const Array &r, Real s) const override
void setLowerBC(const ext::shared_ptr< FdmMesher > &mesher)
FdmSquareRootFwdOp(const ext::shared_ptr< FdmMesher > &mesher, Real kappa, Real theta, Real sigma, Size direction, TransformationType type=Plain)
std::vector< SparseMatrix > toMatrixDecomp() const override
void getCoeffLog(Real &alpha, Real &beta, Real &gamma, Size n) const
ext::shared_ptr< ModTripleBandLinearOp > mapX_
Real upperBoundaryFactor(TransformationType type=Plain) const
void setUpperBC(const ext::shared_ptr< FdmMesher > &mesher)
void setTime(Time t1, Time t2) override
Time is required.
void getCoeffPlain(Real &alpha, Real &beta, Real &gamma, Size n) const
Array apply_mixed(const Array &r) const override
void getCoeffPower(Real &alpha, Real &beta, Real &gamma, Size n) const
Array solve_splitting(Size direction, const Array &r, Real s) const override
const TransformationType transform_
Array apply(const Array &r) const override
void getCoeff(Real &alpha, Real &beta, Real &gamma, Size n) const
Real lowerBoundaryFactor(TransformationType type=Plain) const
Real Time
continuous quantity with 1-year units
std::size_t Size
size of a container
Array Exp(const Array &v)