Docs - requesting peer review

[Luigi Ballabio-4]

Hi all,
        you might have noticed that the QuantLib documentation for release 0.2.0
contains a section which outlines the finite differences framework and
provides a worked example of option pricing---surprisingly enough, it is
titled "The finite differences framework" and it is contained in the "User
manual" part of the documentation.

You all are welcome to read it and send your observations and/or
corrections and/or requests of clarification.

Thanks in advance,
                        Luigi

[Luigi Ballabio-4]

Hi Nicolas,
         I'm sending this to the QuantLib mailing list as well, as it might
be of interest. I do hope you have no problems with that---if you do,
apologies in advance, and blame my rude Italian background.

At 01:58 AM 10/5/01 -0700, you wrote:
>I have a question about Dirichlet boundary
>condition (DBC). The treatement of DBC depends on
>the time scheme, you are using. For forward euler
>(explicit) scheme, it is clear that you set the
>value after each time iteration. For backward
>euler (implicit) scheme, I don't understand what
>is done. Are you also setting the DBC value after
>the time step?

Well, as I didn't write it myself, it took me a few minutes to figure out
what is done---and I did know where to look, too :)

The work is done jointly by three methods, namely,
TridiagonalOperatorCommon::setLowerBC(bc)/setUpperBC(bc) and
TridiagonalOperatorCommon::solveFor(rhs).

The former store the value of the boundary condition and, in the Dirichlet
case, set the upper row of the operator to [1,0,0,....] or the lower row to
[0,0,....,0,1], respectively.

The latter sets the first/last element of the rhs to the stored boundary
condition value(s), then solves the tridiagonal system. This, together with
the operator modification made previously, ensures that the result
satisfies the boundary condition.

Similar tricks are done in the Neumann case.

The drawback I see for the current implementation is that if one creates an
operator with a built-in boundary condition and then sets a Dirichlet or
Neumann one, then it is not possible to revert to the original b.c. unless
one manually sets the first/last row of the operator. With a bit of
housekeeping, this could be done in a future implementation. However, I'm
not at all certain whether the reverting operation should be triggered by
setting a b.c. with type BoundaryCondition::None (meaning "do nothing to
the operator I passed"), or rather some other type we should add to the
enumeration (meaning "the original b.c. the operator had" as opposed to
None meaning "no b.c.").
Feel free to make any suggestions---or even to implement it :)

Hope this helped,
                         Luigi

[Nicolas Di Césaré]

En réponse à Luigi Ballabio <[hidden email]>:

>
> The work is done jointly by three methods, namely,
> TridiagonalOperatorCommon::setLowerBC(bc)/setUpperBC(bc) and
> TridiagonalOperatorCommon::solveFor(rhs).
>
> The former store the value of the boundary condition and, in the
> Dirichlet
> case, set the upper row of the operator to [1,0,0,....] or the lower row
> to
> [0,0,....,0,1], respectively.
>
> The latter sets the first/last element of the rhs to the stored boundary
>
> condition value(s), then solves the tridiagonal system. This, together
> with
> the operator modification made previously, ensures that the result
> satisfies the boundary condition.
>
> Similar tricks are done in the Neumann case.
>
> The drawback I see for the current implementation is that if one creates
> an
> operator with a built-in boundary condition and then sets a Dirichlet or
>
> Neumann one, then it is not possible to revert to the original b.c.
> unless
> one manually sets the first/last row of the operator. With a bit of
> housekeeping, this could be done in a future implementation. However,
> I'm
> not at all certain whether the reverting operation should be triggered
> by
> setting a b.c. with type BoundaryCondition::None (meaning "do nothing to
> the operator I passed"), or rather some other type we should add to the
> enumeration (meaning "the original b.c. the operator had" as opposed to
> None meaning "no b.c.").

> Feel free to make any suggestions---or even to implement it :)

I have worked on a similar design for finite element
method applied to stationnary equations (like Laplace or Poisson equations)
that can be view as backward Euler equation. My conclusion was that boundary
conditions (BC) can modify the equation operator and the right hand side (RHS)
of the equation. Thus, the BC base class interface is :

template <class Vec, class Mat>
class BaseBdyCondition {
public:
  virtual
  void applyTo(const Vec& bdyval, // BC values
               Mat & sparsemat,   // Matrix representation of the operator
               Vec & systemrhs)   // RHS
               = 0;

};

With that interface, I have implemented Neumann, Dirichlet(penalization method)
and Fourrier BC.


The resolution of the problem is made as follow :

  // Define problem type
  Elliptic_PDE Problem(Mesh);

  // Build the equation operator ie set the matrix
  Problem.Equation(laplacian-id);
  // Set the right hand side
  Problem.Set_Rhs(Source);
  // Set Bdy Cond type
  Problem.Bdy_Conditions_Types(Flow_Bdy_Cond);
  // Set Bdy Cond value
  Problem.Bdy_Conditions(Dirichlet_Bdy_Values);
  // Set the resolution method
  Problem.Resolution_Method(gauss);
  // solve the linear system
  Problem.Solve();

This method is wrong for explicit problem where the BC is treated after the
time step u^(k-1) = A.u^k


In conclusion, if we want to allow generic BC to be implemented we have to
   1 - Separate resolution of explicit and implicit schemes.
   2 - use a more generic BC design (like mine ;-) )




--
Nicolas Di Césaré
http://acm.emath.fr/~dicesare

[Luigi Ballabio-4]

Nicolas,

At 03:27 PM 10/5/01 +0200, Nicolas Di C=E9sar=E9 wrote:
>I have worked on a similar design for finite element
>method applied to stationnary equations (like Laplace or Poisson equations)
>that can be view as backward Euler equation. My conclusion was that=
 boundary

>conditions (BC) can modify the equation operator and the right hand side=20
>(RHS)
>of the equation. Thus, the BC base class interface is :
>
>template <class Vec, class Mat>
>class BaseBdyCondition {
>public:
>   virtual
>   void applyTo(const Vec& bdyval, // BC values
>                Mat & sparsemat,   // Matrix representation of the operator
>                Vec & systemrhs)   // RHS
>                =3D 0;
>
>};

I've been thinking along the same lines lately---although without coming to=
=20
your design, which I think covers the whole picture---and I completely=
 agree.

This way we could also pass boundary conditions to the model itself rather=
=20
than setting them to the operator, which kind of makes more sense to me. I=
=20
usually picture to myself boundary conditions as part of the model together=
=20
with the operator, rather than part of the model itself. Of course I'm open=
=20
to corrections on this matter.

Also, this makes me think of passing to the model a function object as well=
=20
which would describe the pay off and would be used to set the initial=20
condition, rather than discretize the latter manually. Now that I think of=
=20
it, that would mean passing the grid as well.

Bottom line: apart from the initial condition thing which needs a bit more=
=20
thinking, if nobody disagrees with your design, I'll put in the todo list.

Bye,
         Luigi

[Luigi Ballabio-4]

Jonathan,
         I'm cc'ing this to the QuantLib users list as well, so that I can
get the proper blame for this blunder.

At 06:42 PM 10/8/01 -0700, Jonathan Knight wrote:
>I have taken a look at the example, and I cannot
>get the interest rate part of the model to work.
>There does not seem to be a way to multiply the
>rate by the identity array in the header files. Am
>  There also seems to be no way to take the
>resulting array and add it to the tridiagonal
>arrays for the rest of the model. Am I missign
>something about how the Identity array works ?

First of all, thanks, This is exactly the reason why I asked for peer
review---I could have looked at that code for hours without seeing this.

You are right: I've been thinking of the Identity as a degenerate
tridiagonal operator and I forgot that it isn't implemented as one. Or
maybe I've been thinking of how it is exported towards Python. Argh!

Well, while I think of a clever way of including the actual Identity class
into the operator algebra---and you are all welcome to contribute
ideas---you can fix the thing as follows:

class TridiagonalIdentity : public TridiagonalOperator {
   public:
     TridiagonalIdentity(int points)
     : TridiagonalOperator(Array(points-1,0.0), // lower diagonal
                           Array(points,1.0),   // diagonal
                           Array(points-1,0.0)) // upper diagonal
     {}
};

class BlackScholesOperator : public TridiagonalOperator {
   as in the docs, replacing Identity<Array>(points) which doesn't work with
   TridiagonalIdentity(points) which does.
};

Apologies,
                 Luigi