Re: [Quantlib-dev] fd questions ctd...

[Peter C. <pca...@vo...>]

Hallo Klaus,

thanks a lot for your answers. I implemented a Neumann condition for the 
multidimensional case. I tested the new class against known solutions 
for the 1d heat equation with Neumann and mixed lower Dirichlet / upper 
Neumann conditions using some of the existing multidim schemes. The 
results suggest that the implementation is ok. Also, the testsuite runs 
without problems which is good I guess.

However, I could not see how to implement the condition by changing the 
final operator directly as done in the 1d framework. Instead I added 
functionality to change the discretization of multidim operators from 
default = FreeBoundary to Neumann by adding a method

void FdmLinearOp::discretization(Size direction, Discretization d);

The enum DiscretizationType in the same class can be used to specify a 
Neumann discretization together with the side (using bitwise or; this 
could be extended for other boundary conditions later). The method calls 
another protected virtual method changeDiscretization(...) which should 
be implemented by user defined operators. I did not do that for the 
existing operators yet (in case I am on the wrong road with my 
approach), but rather added a test operator FdmHeat1dOp (representing 
the pde u_t = \alpha u_xx) which demonstrates the principle. The 
required extension of operators can easily be done being supported by 
extended constructors of the 'basis' operators (first, second and mixed 
second derivatives) allowing for Neumann discretization now. I think I 
did that discretization in the standard way, but it should be cross checked.

I amended the scheme implementations w.r.t. their calls to the 
apply...() methods having the impression that these accounted for the 
Dirichlet case only so far. Needless to say that I am totally not sure 
if I did that correctly. Since I also extended the Dirichlet condition 
to the time (and spatial) dependent case, I inserted the needed 
setTime() calls in the schemes accordingly. In the Dirichlet condition 
itself I implemented the applyBeforeSolving() method because I think it 
is necessary for the time dependent case (and comparisions to known 1d 
heat equation solutions do confirm this).

I removed the FdmBoundaryConditionSet typedef from the multidim 
Dirichlet condition, instead refering to 
OperatorTraits<FdmLinearOp>::bc_set in the schemes, solvers, engines and 
operators now. There is a applyAfterApplying(Real,Real) method in the 
Dirichlet condition (used in some engine I think) which I added to the 
general interface (by default throwing an exception to ensure that only 
the implementation in the Dirichlet boundary is actually used).

Finally I amended the concentrating1dmesher once again because I noted 
that it throws an exception when QL_EXTRA_SAFETY_CHECKS are enabled and 
the concentrating point is equal to one of the endpoints (this having to 
do with an extra check concerning the strict order of x values in 
interpolation.hpp).

I attach my changes (zips with changed files as per directory and new 
files attached directly, all based on ql 1.1) and would be happy if you 
could have a look.

Thank you
Peter


Am 23.04.2012 13:08, schrieb Klaus Spanderen:
> Hi Peter
>
> 1. Thanks for your changes, I've added the new class to the SVN trunk.
>
> 2. At the time being I think the answer is no. We'll have to extend
>   ql/methods/finitedifferences/FiniteDifferenceModel::rollbackImpl
> by the capability to deal with non uniform time grids.
>
> 3. I've used FiniteDifferenceModel::rollback to solve the Fokker-Planck
> equation. Even though the name of the method is "rollback", if "from is
> smaller than to" the methods rolls forward. You'll have to remove line 95 and
> 96
>
>              QL_REQUIRE(from>= to,
>                         "trying to roll back from "<<  from<<  " to "<<  to);
>
> from finitedifferencemodel.hpp to get it working.
>
> 4. The multidim frameworks supports "free boundary condition" (default if no
> BC is given) and Dirichlet BC. Your are right, the Neumann condition was
> never implemented. The boundary conditions are hardcoded in some classes
> because otherwise users might supply the one dimensional classes NeumannBC or
> DirichletBC, which don't work in the multidimensional case. If you have a
> multidimensional version of NeumannBC running we can change the interfaces.
>
> regards
>   Klaus
>
> On Sunday 22 April 2012 19:27:12 Peter Caspers wrote:
>> Hi,
>>
>> may I ask one more question please:
>>
>> 4. it seems there is only a Dirichlet boundary condition implemented in
>> the multi dimensional context. I want a Neumann condition and I think I
>> managed to implement a version, but the boundary conditions are
>> hardcoded as FdmDirichletBoundary in many classes, even the typedef for
>> the FdmBoundaryConditionSet is vector<shared_ptr<FdmDirichletBoundary>>.
>> Am I missing something here? How else would I be supposed to add new
>> boundary conditions?
>>
>> As for 1. below I extended the concentrating mesher by a flag that
>> allows the central point to be forced into the mesh (using a piecewise
>> linear transformation of the generating uniform grid, see e.g. Iain
>> Clarke, FX Option Pricing, ch. ...). To set up non uniform grids with
>> more than one concentrating point I added the glued1dmesher. If you
>> consider these contributions useful, please add them to the library.
>>
>> 2 and 3 have obvious workarounds (rewriting the pde as backward and
>> doing the time steps manually one by one). Still the implementation of
>> forward operators becomes less readable by this implicit (yet trivial)
>> transformation and possibly it may be useful in general to have a non
>> uniform time grid in the solvers, e.g. via a transformation [0,1] ->
>> [0,1], u ->  pow(u,alpha), alpha>  0 and a linear one [0,1] ->  [0,T] (cf.
>> same reference as above). If considered useful, I'd be happy to do these
>> extensions to the FdmBackwardSolver class.
>>
>> Regards
>> Peter
>>
>> -------- Original-Nachricht --------
>> Betreff: 	[Quantlib-users] fd questions
>> Datum: 	Thu, 12 Apr 2012 13:37:25 +0200
>> Von: 	Peter Caspers<pca...@vo...>
>> An: 	qua...@li...
>>
>>
>>
>> Hi,
>>
>> I just started to use the ql 1.1 / finitedifferences framework and have
>> a couple of (probably very basic) questions:
>>
>> 1. Is there a way to specify a or even several mandatory point(s) in the
>> meshers?
>> 2. Can I use a non uniform time grid in the solver?
>> 3. Is there a forward solver ?
>>
>> Thank you
>> Peter
>>
>>
>>
>>
>> ---------------------------------------------------------------------------
>> --- For Developers, A Lot Can Happen In A Second.
>> Boundary is the first to Know...and Tell You.
>> Monitor Your Applications in Ultra-Fine Resolution. Try it FREE!
>> http://p.sf.net/sfu/Boundary-d2dvs2
>> _______________________________________________
>> QuantLib-users mailing list
>> Qua...@li...
>> https://lists.sourceforge.net/lists/listinfo/quantlib-users
>
> ------------------------------------------------------------------------------
> For Developers, A Lot Can Happen In A Second.
> Boundary is the first to Know...and Tell You.
> Monitor Your Applications in Ultra-Fine Resolution. Try it FREE!
> http://p.sf.net/sfu/Boundary-d2dvs2
> _______________________________________________
> QuantLib-dev mailing list
> Qua...@li...
> https://lists.sourceforge.net/lists/listinfo/quantlib-dev