Re: [Quantlib-dev] [QuantLib-svn] SF.net SVN: quantlib:[18270] trunk/QuantLib

SourceForge Headquarters 1320 Columbia Street Suite 310 San Diego, CA 92101 +1 (858) 422-6466

Last one :)

Is the order of convergence (or a lack of one) known when an
extrapolation is created, given h and f?  If so, we could pass it to
the constructor and leave the interface for calling the formula just
as operator()(t) so that we don't have to pass the order again and
again,

We might also overload it with operator()(t,s) if we didn't pass the
order of convergence.  Just a thought, no big deal anyway---it's just
that formulaUnknownConvergence looks awfully long...

Luigi

On Sun, Jun 10, 2012 at 10:56 AM,  <kla...@us...> wrote:
> Revision: 18270
>          http://quantlib.svn.sourceforge.net/quantlib/?rev=18270&view=rev
> Author:   klausspanderen
> Date:     2012-06-10 08:56:09 +0000 (Sun, 10 Jun 2012)
> Log Message:
> -----------
> added Richardson extrapolation
>
> Modified Paths:
> --------------
>    trunk/QuantLib/ql/math/Makefile.am
>    trunk/QuantLib/ql/math/all.hpp
>    trunk/QuantLib/test-suite/interpolations.cpp
>    trunk/QuantLib/test-suite/interpolations.hpp
>
> Added Paths:
> -----------
>    trunk/QuantLib/ql/math/richardsonextrapolation.cpp
>    trunk/QuantLib/ql/math/richardsonextrapolation.hpp
>
> Added: trunk/QuantLib/ql/math/richardsonextrapolation.cpp
> ===================================================================
> --- trunk/QuantLib/ql/math/richardsonextrapolation.cpp                          (rev 0)
> +++ trunk/QuantLib/ql/math/richardsonextrapolation.cpp  2012-06-10 08:56:09 UTC (rev 18270)
> @@ -0,0 +1,73 @@
> +/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
> +
> +/*
> + Copyright (C) 2012 Klaus Spanderen
> +
> + This file is part of QuantLib, a free-software/open-source library
> + for financial quantitative analysts and developers - http://quantlib.org/
> +
> + QuantLib is free software: you can redistribute it and/or modify it
> + under the terms of the QuantLib license.  You should have received a
> + copy of the license along with this program; if not, please email
> + <qua...@li...>. The license is also available online at
> + <http://quantlib.org/license.shtml>.
> +
> + This program is distributed in the hope that it will be useful, but WITHOUT
> + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
> + FOR A PARTICULAR PURPOSE.  See the license for more details.
> +*/
> +
> +/*! \file richardsonextrapolation.cpp
> +*/
> +
> +#include <ql/errors.hpp>
> +#include <ql/math/solvers1d/brent.hpp>
> +#include <ql/math/richardsonextrapolation.hpp>
> +
> +#include <cmath>
> +
> +namespace QuantLib {
> +    namespace {
> +        class RichardsonEqn {
> +          public:
> +            RichardsonEqn(Real fh, Real ft, Real fs, Real t, Real s)
> +            : fh_(fh), ft_(ft), fs_(fs), t_(t), s_(s) { }
> +
> +            Real operator()(Real k) const {
> +                return      ft_ + (ft_-fh_)/(std::pow(t_, k)-1.0)
> +                        - ( fs_ + (fs_-fh_)/(std::pow(s_, k)-1.0));
> +            }
> +          private:
> +            const Real fh_, ft_, fs_, t_, s_;
> +        };
> +
> +    }
> +
> +    RichardsonExtrapolation::RichardsonExtrapolation(
> +        Real h, const boost::function<Real (Real)>& f)
> +    : h_(h), fh_(f(h)), f_(f) { }
> +
> +
> +    Real RichardsonExtrapolation::formula(Real t, Real k) const {
> +
> +        QL_REQUIRE(t > 1, "scaling factor must be greater than 1");
> +        const Real tk = std::pow(t,k);
> +
> +        return (tk*f_(h_/t)-fh_)/(tk-1.0);
> +    }
> +
> +    Real RichardsonExtrapolation::formulaUnknownConvergence(Real t, Real s)
> +    const {
> +        QL_REQUIRE(t > 1 && s > 1, "scaling factors must be greater than 1");
> +
> +        const Real ft = f_(h_/t);
> +        const Real fs = f_(h_/s);
> +
> +        const Real k = Brent().solve(RichardsonEqn(fh_, ft, fs, t, s),
> +                                     1e-8, 0.05, 10);
> +
> +        const Real ts = std::pow(s,k);
> +
> +        return (ts*fs-fh_)/(ts-1.0);
> +    }
> +}
>
>
> Property changes on: trunk/QuantLib/ql/math/richardsonextrapolation.cpp
> ___________________________________________________________________
> Added: svn:mime-type
>   + text/plain
>
> Added: trunk/QuantLib/ql/math/richardsonextrapolation.hpp
> ===================================================================
> --- trunk/QuantLib/ql/math/richardsonextrapolation.hpp                          (rev 0)
> +++ trunk/QuantLib/ql/math/richardsonextrapolation.hpp  2012-06-10 08:56:09 UTC (rev 18270)
> @@ -0,0 +1,62 @@
> +/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
> +
> +/*
> + Copyright (C) 2012 Klaus Spanderen
> +
> + This file is part of QuantLib, a free-software/open-source library
> + for financial quantitative analysts and developers - http://quantlib.org/
> +
> + QuantLib is free software: you can redistribute it and/or modify it
> + under the terms of the QuantLib license.  You should have received a
> + copy of the license along with this program; if not, please email
> + <qua...@li...>. The license is also available online at
> + <http://quantlib.org/license.shtml>.
> +
> + This program is distributed in the hope that it will be useful, but WITHOUT
> + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
> + FOR A PARTICULAR PURPOSE.  See the license for more details.
> +*/
> +
> +/*! \file richardsonextrapolation.hpp
> +*/
> +
> +#ifndef quantlib_richardson_extrapolation_hpp
> +#define quantlib_richardson_extrapolation_hpp
> +
> +#include <ql/types.hpp>
> +
> +#include <boost/function.hpp>
> +
> +namespace QuantLib {
> +
> +    //! Richardson Extrapolatio
> +    /*! References:
> +
> +        http://en.wikipedia.org/wiki/Richardson_extrapolation
> +     */
> +
> +    class RichardsonExtrapolation {
> +      public:
> +        RichardsonExtrapolation(Real h,
> +                                const boost::function<Real (Real)>& f);
> +
> +        /*! Extrapolation for known order of convergence
> +            \param  step size for extrapolation. Is often 2
> +            \param k order of convergence
> +        */
> +        Real formula(Real t, Real k) const;
> +
> +        /*! Extrapolation for unknown order of convergence
> +            \param t first step size for the extrapolation.
> +            \param s second step size for the extrapolation
> +        */
> +        Real formulaUnknownConvergence(Real t, Real s) const;
> +
> +      private:
> +        const Real h_;
> +        const Real fh_;
> +        const boost::function<Real (Real)> f_;
> +    };
> +}
> +
> +#endif
>
>
> Property changes on: trunk/QuantLib/ql/math/richardsonextrapolation.hpp
> ___________________________________________________________________
> Added: svn:mime-type
>   + text/plain
>
> Modified: trunk/QuantLib/test-suite/interpolations.cpp
> ===================================================================
> --- trunk/QuantLib/test-suite/interpolations.cpp        2012-05-30 22:32:14 UTC (rev 18269)
> +++ trunk/QuantLib/test-suite/interpolations.cpp        2012-06-10 08:56:09 UTC (rev 18270)
> @@ -36,6 +36,7 @@
>  #include <ql/math/integrals/simpsonintegral.hpp>
>  #include <ql/math/kernelfunctions.hpp>
>  #include <ql/math/functional.hpp>
> +#include <ql/math/richardsonextrapolation.hpp>
>  #include <ql/math/randomnumbers/sobolrsg.hpp>
>  #include <ql/math/optimization/levenbergmarquardt.hpp>
>  #include <boost/assign/std/vector.hpp>
> @@ -1610,6 +1611,35 @@
>     }
>  }
>
> +namespace {
> +    Real f(Real h) {
> +        return std::pow( 1.0 + h, 1/h);
> +    }
> +}
> +
> +void InterpolationTest::testRichardsonExtrapolation() {
> +    BOOST_MESSAGE("Testing Richardson Extrapolation...");
> +
> +    /* example taken from
> +     * http://www.ipvs.uni-stuttgart.de/abteilungen/bv/lehre/
> +     *      lehrveranstaltungen/vorlesungen/WS0910/
> +     *      NSG_termine/dateien/Richardson.pdf
> +     */
> +
> +    const RichardsonExtrapolation extrap(0.1, f);
> +
> +    const Real stepSize = 2.0;
> +    const Real orderOfConvergence = 1.0;
> +
> +    const Real tol = 0.00001;
> +    const Real expected = 2.71285;
> +    const Real calculated = extrap.formula(stepSize, orderOfConvergence);
> +
> +    if (std::fabs(expected-calculated) > tol) {
> +        BOOST_ERROR("failed to reproduce Richardson extrapolation");
> +    }
> +}
> +
>  test_suite* InterpolationTest::suite() {
>     test_suite* suite = BOOST_TEST_SUITE("Interpolation tests");
>
> @@ -1636,6 +1666,8 @@
>     suite->add(QUANTLIB_TEST_CASE(
>                               &InterpolationTest::testKernelInterpolation2D));
>     suite->add(QUANTLIB_TEST_CASE(&InterpolationTest::testBicubicDerivatives));
> +    suite->add(QUANTLIB_TEST_CASE(
> +                            &InterpolationTest::testRichardsonExtrapolation));
>
>     return suite;
>  }
>
> Modified: trunk/QuantLib/test-suite/interpolations.hpp
> ===================================================================
> --- trunk/QuantLib/test-suite/interpolations.hpp        2012-05-30 22:32:14 UTC (rev 18269)
> +++ trunk/QuantLib/test-suite/interpolations.hpp        2012-06-10 08:56:09 UTC (rev 18270)
> @@ -46,6 +46,8 @@
>     static void testKernelInterpolation();
>     static void testKernelInterpolation2D();
>     static void testBicubicDerivatives();
> +    static void testRichardsonExtrapolation();
> +
>     static boost::unit_test_framework::test_suite* suite();
>  };
>
>
> This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
>
>
> ------------------------------------------------------------------------------
> Live Security Virtual Conference
> Exclusive live event will cover all the ways today's security and
> threat landscape has changed and how IT managers can respond. Discussions
> will include endpoint security, mobile security and the latest in malware
> threats. http://www.accelacomm.com/jaw/sfrnl04242012/114/50122263/
> _______________________________________________
> QuantLib-cvs mailing list
> Qua...@li...
> https://lists.sourceforge.net/lists/listinfo/quantlib-cvs