Login  Register

RE: Docs - requesting peer review

Posted by Luigi Ballabio-4 on Oct 05, 2001; 4:15am
URL: http://quantlib.414.s1.nabble.com/Docs-requesting-peer-review-tp1768p1770.html

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