mapcubi.h

Go to the documentation of this file.

/// @addtogroup unifexp
/// @{

/////////////////////////////////////////////////////////////////////////////
///
/// @file mapcubi.h
///
/// This file contains the definition of the cubic map -x^3 + 3x - a
/// on [-4,4] with the usage of interval arithmetic
/// for rigorous computations.
///
/// @author Pawel Pilarczyk
///
/////////////////////////////////////////////////////////////////////////////

// Copyright (C) 2007 by Pawel Pilarczyk.
//
// This file is part of my research program package.  This is free software;
// you can redistribute it and/or modify it under the terms of the GNU
// General Public License as published by the Free Software Foundation;
// either version 2 of the License, or (at your option) any later version.
//
// This software 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
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License along
// with this software; see the file "license.txt".  If not, write to the
// Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
// MA 02111-1307, USA.

// Started on August 24, 2007. Last revision: February 3, 2013.

#ifndef _mapcubi_h_
#define _mapcubi_h_

#include <string>
#include <iostream>
#include "maptype.h"
#include "intvcapd.h"


namespace unifexp {

// --------------------------------------------------
// ----------------- the cubic map ------------------
// --------------------------------------------------

/// This class defines the cubic map -x^3 + 3x - a on [-4,4]
/// with the usage of interval arithmetic.
/// It is suitable for rigorous computations.
/// Good values of the parameter are between 0 and 2,
/// more or less (the optimal value is 2, I think).
template <class numType>
class mapCubicIntv: public mapType<numType>
{
public:
        /// The constructor.
        mapCubicIntv ();

        /// Returns the name of the object.
        std::string name () const;

        /// Returns the number of critical points.
        int countCritical () const;

        /// Returns the subsequent critical points.
        numType criticalPoint (int n) const;

        /// Returns the left bound of the domain of the map.
        numType leftBound () const;

        /// Returns the right bound of the domain of the map.
        numType rightBound () const;

        /// Computes an enclosure of the image of the given interval.
        void image (const numType &x1, const numType &x2,
                numType &y1, numType &y2) const;

        /// Computes the minimal log of the derivative over those points
        /// in the interval [x1,x2] whose images may fall into [y1,y2]
        numType minLogDerivative (const numType &x1, const numType &x2,
                const numType &y1, const numType &y2) const;

}; /* mapCubicIntv */

// --------------------------------------------------

template <class numType>
inline mapCubicIntv<numType>::mapCubicIntv ()
{
        return;
} /* mapCubicIntv::mapCubicIntv */

template <class numType>
inline std::string mapCubicIntv<numType>::name () const
{
        return std::string ("cubic-intv");
} /* mapCubicIntv::name */

template <class numType>
inline int mapCubicIntv<numType>::countCritical () const
{
        return 2;
} /* mapCubicIntv::countCritical */

template <class numType>
inline numType mapCubicIntv<numType>::criticalPoint (int n) const
{
        return n ? 1 : -1;
} /* mapCubicIntv::criticalPoint */

template <class numType>
inline numType mapCubicIntv<numType>::leftBound () const
{
        return -4;
} /* mapCubicIntv::leftBound */

template <class numType>
inline numType mapCubicIntv<numType>::rightBound () const
{
        return 4;
} /* mapCubicIntv::rightBound */

template <class numType>
inline void mapCubicIntv<numType>::image (const numType &x1,
        const numType &x2, numType &y1, numType &y2) const
{
        if (x2 < x1)
                throw "Image computation: Wrong interval for 'x'.";
        IntervalType x1i (x1, x1);
        IntervalType x2i (x2, x2);
        if ((x2 <= -1) || (1 <= x1))
        {
                y1 = (-(x2i * x2i - 3) * x2i - this -> paramMax).
                        leftBound ();
                y2 = (-(x1i * x1i - 3) * x1i - this -> paramMin).
                        rightBound ();
        }
        else if ((-1 <= x1) && (x2 <= 1))
        {
                y1 = (-(x1i * x1i - 3) * x1i - this -> paramMax).
                        leftBound ();
                y2 = (-(x2i * x2i - 3) * x2i - this -> paramMin).
                        rightBound ();
        }
        else if (x2 <= 1)
        {
                const IntervalType y3 = -(x1i * x1i - 3) * x1i;
                const IntervalType y4 = -(x2i * x2i - 3) * x2i;
                y1 = (IntervalType (-2, -2) - this -> paramMax).
                        leftBound ();
                y2 = (((y3. rightBound () < y4. rightBound ()) ? y4 : y3) -
                        this -> paramMin). rightBound ();
        }
        else if (-1 <= x1)
        {
                const IntervalType y3 = -(x1i * x1i - 3) * x1i;
                const IntervalType y4 = -(x2i * x2i - 3) * x2i;
                y1 = (((y3. leftBound () < y4. leftBound ()) ? y3 : y4) -
                        this -> paramMax). leftBound ();
                y2 = (IntervalType (4, 4) - this -> paramMin). rightBound ();
        }
        else
                throw "Cubic map: Too large domain interval.";
        resetRounding ();
        return;
} /* mapCubicIntv::image */

template <class numType>
inline numType mapCubicIntv<numType>::minLogDerivative (const numType &x1,
        const numType &x2, const numType &y1, const numType &y2) const
{
        // make sure the input data is correct
        if (x2 < x1)
                throw "MinLogDerivative: Wrong interval for 'x'.";
        if (y2 < y1)
                throw "MinLogDerivative: Wrong interval for 'y'.";

        // NOTE: Since it is not so easy to compute the preimage of the
        // interval [y1,y2] with respect to the cubic polynomial map,
        // the lower bound for the log of the derivative is computed
        // solely based on the interval [x1,x2]
        IntervalType x1i (x1, x1);
        IntervalType x2i (x2, x2);
        numType result;
        if (x2 < -1)
        {
                result = log (3 * x2i * x2i + 3). leftBound ();
        }
        else if (1 < x1)
        {
                result = log (3 * x1i * x1i + 3). leftBound ();
        }
        else if ((-1 < x1) && (x2 < 1))
        {
                if (0 < x1)
                        result = log (-3 * x2i * x2i + 3). leftBound ();
                else if (x2 < 0)
                        result = log (-3 * x1i * x1i + 3). leftBound ();
                else
                {
                        result = log (-3 * ((x1 < x2) ? (x2i * x2i) :
                                (x1i * x1i)) + 3). leftBound ();
                }
        }
        else
        {
                throw "Log of interval containing zero.";
                result = 0;
        }
        resetRounding ();
        return result;
} /* mapCubicIntv::minLogDerivative */


} // namespace unifexp

#endif // _mapcubi_h_

/// @}