Hello everybody.

This week’s content is another installment in the series on finite-difference methods.

Subscribe to my Substack to receive my posts in your inbox, or follow me on Twitter or LinkedIn if you want to be notified of new posts, or subscribe via RSS if you’re the tech type: the buttons for all that are in the footer. Also, I’m available for training, both online and (when possible) on-site: visit my Training page for more information.

Boundary conditions

As I said, there’s more going on in the step method. To begin with, we’ll have to enforce boundary conditions at either end of the grid; for instance, we might want to ensure that an option has price close to 0 at the grid boundary which is deep out of the money and derivative close to 1 at the boundary which is deep in the money. Such conditions are modeled by the BoundaryCondition class, shown in the listing below.

    template <class Operator>
    class BoundaryCondition {
      public:
        typedef Operator operator_type;
        typedef typename Operator::array_type array_type;

        enum Side { None, Upper, Lower };

        virtual ~BoundaryCondition() {}

        virtual void applyBeforeApplying(operator_type&) const = 0;
        virtual void applyAfterApplying(array_type&) const = 0;
        virtual void applyBeforeSolving(operator_type&,
                                        array_type& rhs) const = 0;
        virtual void applyAfterSolving(array_type&) const = 0;
    };

As you see, the class uses runtime polymorphism. This is because we wanted to store collections of them, possibly of different types (you might remember the bc_set typedef from the listing of the MixedScheme class), and most suitable containers are homogeneous and thus require a common base class (std::tuple wasn’t even a gleam in the standard committee’s eye, and handling it requires template metaprogramming, which I’d try to avoid unless necessary—even in this time and age) We could have used std::pair if we had decided to limit ourselves to the 1-D case; we didn’t, though, and in hindsight I think it avoided adding more complexity to the framework.

The interface of the class is, well, peculiar. Glossing over the Side enumeration (which is meant to specify a grid side in 1-D models, and thus too specific), methods like applyBeforeApplying are plausible contenders for the worst name in the library. The idea here is that the boundary condition can be enforced, or “applied”, in different ways. One can wait for the operator’s applyTo or solveFor to execute, and then modify the result; this would be implemented by applyAfterApplying and applyAfterSolving, respectively. Another possibility is to set up the operator and the input array beforehand, so that the result of the operation will satisfy the condition; this is done in applyBeforeApplying and applyBeforeSolving. (For some reason, the former only takes the operator. It might have been an oversight, due to the fact that the existing conditions didn’t require to modify the array.)

As an example, look at the DirichletBC class in the listing below. It implements a simple Dirichlet boundary condition, i.e., one in which the function value at a given end of the grid must equal a given constant value.

    class DirichletBC
        : public BoundaryCondition<TridiagonalOperator> {
      public:
        DirichletBC(Real value, Side side);

        void applyBeforeApplying(TridiagonalOperator& L) const {
            switch (side_) {
              case Lower:
                L.setFirstRow(1.0,0.0);
                break;
              case Upper:
                L.setLastRow(0.0,1.0);
                break;
            }
        }
        void applyAfterApplying(Array& u) const {
            switch (side_) {
              case Lower:
                u[0] = value_;
                break;
              case Upper:
                u[u.size()-1] = value_;
                break;
            }
        }
        void applyBeforeSolving(TridiagonalOperator& L,
                                Array& rhs) const {
            switch (side_) {
              case Lower:
                L.setFirstRow(1.0,0.0);
                rhs[0] = value_;
                break;
              case Upper:
                L.setLastRow(0.0,1.0);
                rhs[rhs.size()-1] = value_;
                break;
            }
        }
        void applyAfterSolving(Array&) const {}
    };

As you can see, it is no longer a class template; when implemented, boundary conditions are specialized for a given operator—not surprisingly, because they have to know how to access and modify it. In this case, we’re using the tridiagonal operator I described in an earlier subsection.

The constructor is simple enough: it takes the constant value and the side of the grid to which it must be applied, and stores them.

As I said, the applyBeforeApplying method must set up the operator L so that the result of L.apply(u), or \( \mathbf{u^{\prime}} = L \cdot \mathbf{u} \), satisfies the boundary condition. To do this, it sets the first (last) row to the corresponding row of the identity matrix; this ensures that the first (last) element of \( \mathbf{u^{\prime}} \) equal the corresponding element of \( \mathbf{u} \). Given that the array satisfied the boundary condition at the previous step, it keeps satisfying it at this step. Relying on this is a bit unsafe (for instance, it would break if the value changed at different times, or if any other event modified the array values) but it’s the best one can do without accessing the input array. The applyAfterApplying method is simpler and safer: it just sets the value of the output array directly.

The applyBeforeSolving method works as its sibling applyBeforeApplying, but it can also access the input array, so it can ensure that it contains the correct value. Finally, and surprisingly enough, as of this writing the applyAfterSolving method is empty. When we wrote this class, we probably relied on the fact that applyBeforeSolving would also be called, which doesn’t strike me as a very good idea in hindsight: we’re not sure of what methods any given evolution scheme may call (this also makes applyBeforeApplying less safe than it seems).

The reason this worked is the way the MixedScheme::step method is implemented, as seen in the next listing. At each step, all boundary conditions are enforced first before and after applying the explicit part of the scheme, and then before and after solving for the implicit part. In the case of the Dirichlet condition, this causes applyAfterApplying to fix any problem that applyBeforeApplying might have introduced and also ensures that applyBeforeSolving does the work even if applyAfterSolving doesn’t.

If a boundary condition class implemented all of its methods correctly, this would be redundant; enforcing the condition just after the operations would be enough. However, we probably assumed that some conditions couldn’t (or just didn’t) implement all of them, and went for the safe option.

    template <class Operator>
    void MixedScheme<Operator>::step(array_type& a, Time t) {
        if (theta_ != 1.0) { // there is an explicit part
            for (Size i=0; i<bcs_.size(); i++)
                bcs_[i]->applyBeforeApplying(explicitPart_);
            a = explicitPart_.applyTo(a);
            for (Size i=0; i<bcs_.size(); i++)
                bcs_[i]->applyAfterApplying(a);
        }
        if (theta_ != 0.0) { // there is an implicit part
            for (Size i=0; i<bcs_.size(); i++)
                bcs_[i]->applyBeforeSolving(implicitPart_, a);
            implicitPart_.solveFor(a, a);
            for (Size i=0; i<bcs_.size(); i++)
                bcs_[i]->applyAfterSolving(a);
        }
    }

If we ever set out to revamp this part of the library (not that likely, given the new framework) this is something we should get a look at. The problem to solve, unfortunately, is not an easy one: to call only the methods that are required by the scheme and supported by the boundary condition, we’d need to know the specific type of both.

Inverting the dependency and passing the scheme to the boundary condition, instead of the other way around, wouldn’t work: the boundary condition wouldn’t know if the scheme is calling applyTo or solveFor, and wouldn’t be able to act inside the scheme’s step method (as in MixedScheme::step, where the boundary conditions are enforced between the implicit and explicit parts of the step.

A mediator or some kind of visitor might work, but the complexity of the code would increase (especially if the number of conditions grow). Had they been available in C++, multimethods might have been the best choice; but even those might not be usable for multiple conditions. All in all, the current workaround of calling all boundary condition method might be redundant but at least works correctly. The only alternative I can think of would be to require all conditions to be enforced after the operation, remove the applyBefore... methods and be done with it.

Enough with boundary conditions for now. With the pieces we have so far, we can take a payoff at \( t=T \) and evolve it back step by step to \( t=0 \)—if nothing happens in between, that is. But that’s not often the case, which is my cue for the next post.