Boundary conditions (BC) are an essential part of a finite-difference model.
In the current implementation, BC:s are implemented as a simple structure with a enumerated type tag and a value. The logic of the boundary condition itself is implemented in the TridiagonalOperator class. This has drawbacks in terms of both model logic and extensibility.
A new implementation is proposed which encapsulates the logic of the BC into a polymorphic BoundaryCondition class which is part of the finite-difference model as much as the differential operator and the time scheme.
The present implementation of FiniteDifferenceMethod contains an Evolver, which in turn contains an Operator. BC:s are set to the operator before the latter is passed to the model which in turn passes it to the evolver. Setting a BC immediately modifies an operator; also, BC:s are stored in the operator and applied to the input and/or output array when the applyTo() or solveFor() methods are invoked.
There is no clean encapsulation of the BC logic. The operator and the input/output arrays are modified at different times and different places in the code. Also, there is no clear separation between the BC logic and the linear algebra calculations which form the core of the applyTo() and solveFor() methods.
The code which applies the BC:s is a switch construct which executes different code depending on the enumerated BC type tag. This seriously compromises the extensibility of the code, e.g., the addition of new BC types to the framework.
Another logic flaw is present, namely, BC:s are set and modifications are applied to the operator D which is passed to the model. Instead, they should only be applied to the operator which transforms the array at a given time step into the array at the next step. Such operator is not D itself but rather I+dt*D or a similar one depending on the chosen time scheme. For this same reason, care must be taken that the implementation of the operator algebra propagates the BC:s from D to (dt*D) to (I+(dt*D)).
A polymorphic BoundaryCondition class is defined which encapsulates the BC logic. Its interface is defined so that it should be able to manage different kind of required modifications, namely, ones that should be applied before and/or after the application of the operator or its inverse. Such methods act simultaneously on the operator and the right-hand side of the equation. Also, a setTime() interface is provided which allows to manage time-dependent BC:s.
The BC:s are passed to the model together with the operator and are stored into the evolver. It is the responsibility of the evolver to call the appropriate BC methods inside its step() method.
Finally, the TridiagonalOperator would not be concerned with BC:s. Its applyTo() and solveFor() methods would only implement the corresponding linear algebra calculations.
The framework would be implemented as in the following pseudo-code:
class FiniteDifferenceMethod {
public:
FiniteDifferenceMethod(const Operator& L, ??? boundaryConditions)
: evolver_(L,boundaryConditions) { ... }
...
private:
Evolver evolver_;
};
class ExplicitEuler {
public:
ExplicitEuler(const Operator& L, ??? boundaryConditions)
: L_(L), boundaryConditions_(boundaryConditions) { ... }
void step(Array& a, Time t) {
... // call L.setTime(t) if needed
// actual operator is O = I-dt*L
for every boundary condition
bc.setTime(t);
bc.applyBeforeApplying(O);
a = O.applyTo(a);
for every boundary condition
bc.applyAfterApplying(a);
}
};
class ImplicitEuler {
public:
ImplicitEuler(const Operator& L, ??? boundaryConditions)
: L_(L), boundaryConditions_(boundaryConditions) { ... }
void step(Array& a, Time t) {
... // call L.setTime(t-dt) if needed
// actual operator is O = I+dt*L
for every boundary condition
bc.setTime(t-dt);
bc.applyBeforeSolving(O,a);
a = O.solveFor(a);
for every boundary condition
bc.applyAfterSolving(a);
}
};
class CrankNicolson {
public:
CrankNicolson(const Operator& L, ??? boundaryConditions)
: L_(L), boundaryConditions_(boundaryConditions) { ... }
void step(Array& a, Time t) {
... // call L.setTime(t) if needed
// actual operators are O1 = I-(dt/2)*L and O2 = I+(dt/2)*L
for every boundary condition
bc.setTime(t);
bc.applyBeforeApplying(O1);
a = O1.applyTo(a);
for every boundary condition
bc.applyAfterApplying(a);
for every boundary condition
bc.setTime(t-dt);
bc.applyBeforeSolving(O2,a);
a = O2.solveFor(a);
for every boundary condition
bc.applyAfterSolving(a);
}
};
The missing piece in the above code concerns the way to pass the BC:s to a model. This must take into account that we do not know in advance their number even in 1-D models, much less in 2-D or higher.
One possibility would be to pass a container of BC:s, e.g., a std::vector or a std::list. In this case, the constructor of a model would be implemented as
class FiniteDifferenceMethod {
public:
FiniteDifferenceMethod(const Operator& L,
const std::list<Handle<BoundaryCondition> >& bc)
: evolver_(L,bc) { ... }
...
};
Another would be to pass to the model one BC which represents the set of all BC:s which are local to a part of the solution. The Composite pattern could be used to define such set as
class CompositeBC : public BoundaryCondition {
public:
CompositeBC(const std::list<Handle<BoundaryCondition> >& bc =
std::list<Handle<BoundaryCondition> >())
: bc_(bc) {}
void add(const Handle<BoundaryCondition>& c) {
bc_.push_back(c);
}
void setTime(Time t) {
for (every BC in the list)
i->setTime(t);
}
void applyBeforeApplying(Operator& O) {
for (every BC in the list)
i->applyBeforeApplying(O);
}
...
};
In this case, the constructor of a model would be implemented as:
class FiniteDifferenceMethod {
public:
FiniteDifferenceMethod(const Operator& L,
const Handle<BoundaryCondition>& bc)
: evolver_(L,bc) { ... }
...
};
Many thanks to Nicolas Di Césaré for the fruitful discussion---which is in fact still ongoing and still fruitful.
Feedback on the above proposal should be posted to the QuantLib-dev mailing list.