mapunimi.h

Go to the documentation of this file.
00001 /// @addtogroup unifexp
00002 /// @{
00003 
00004 /////////////////////////////////////////////////////////////////////////////
00005 ///
00006 /// @file mapunimi.h
00007 ///
00008 /// This file contains the definition of the unimodal map "2|x|^gamma - a"
00009 /// on [-1,1] with the usage of interval arithmetic
00010 /// for rigorous computations.
00011 ///
00012 /// @author Pawel Pilarczyk
00013 ///
00014 /////////////////////////////////////////////////////////////////////////////
00015 
00016 // Copyright (C) 2007 by Pawel Pilarczyk.
00017 //
00018 // This file is part of my research program package.  This is free software;
00019 // you can redistribute it and/or modify it under the terms of the GNU
00020 // General Public License as published by the Free Software Foundation;
00021 // either version 2 of the License, or (at your option) any later version.
00022 //
00023 // This software is distributed in the hope that it will be useful,
00024 // but WITHOUT ANY WARRANTY; without even the implied warranty of
00025 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00026 // GNU General Public License for more details.
00027 //
00028 // You should have received a copy of the GNU General Public License along
00029 // with this software; see the file "license.txt".  If not, write to the
00030 // Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
00031 // MA 02111-1307, USA.
00032 
00033 // Started on August 23, 2007. Last revision: August 24, 2007.
00034 
00035 #ifndef _mapunimi_h_
00036 #define _mapunimi_h_
00037 
00038 #include <string>
00039 #include <iostream>
00040 #include <sstream>
00041 #include "maptype.h"
00042 
00043 // the "interval" package from the CAPD library
00044 #include "interval/doubleInterval.h"
00045 
00046 
00047 namespace unifexp {
00048 
00049 // --------------------------------------------------
00050 // ---------------- the unimodal map ----------------
00051 // --------------------------------------------------
00052 
00053 /// This class defines the unimodal map "2 |x|^gamma - a" on [-1,1]
00054 /// with the usage of interval arithmetic.
00055 /// It is suitable for rigorous computations.
00056 /// Recommended values of gamma are between 1 and 3.
00057 /// Valid values of the parameter are between 0.5+ and 1 (optimal value: 1).
00058 template <class numType>
00059 class mapUnimodalIntv: public mapType<numType>
00060 {
00061 public:
00062         /// The constructor. Provide the numerator and the denominator
00063         /// for the exponent gamma to make sure its value is computed
00064         /// in a rigorous way as an interval.
00065         mapUnimodalIntv (const numType &_gamma1, const numType &_gamma2 = 1);
00066 
00067         /// Returns the name of the object.
00068         std::string name () const;
00069 
00070         /// Returns the number of critical points.
00071         int countCritical () const;
00072 
00073         /// Returns the subsequent critical points.
00074         numType criticalPoint (int n) const;
00075 
00076         /// Returns the left bound of the domain of the map.
00077         numType leftBound () const;
00078 
00079         /// Returns the right bound of the domain of the map.
00080         numType rightBound () const;
00081 
00082         /// Computes an enclosure of the image of the given interval.
00083         void image (const numType &x1, const numType &x2,
00084                 numType &y1, numType &y2) const;
00085 
00086         /// Computes the minimal log of the derivative over those points
00087         /// in the interval [x1,x2] whose images may fall into [y1,y2]
00088         numType minLogDerivative (const numType &x1, const numType &x2,
00089                 const numType &y1, const numType &y2) const;
00090 
00091 private:
00092         /// The exponent of the map.
00093         interval gamma;
00094 
00095         /// The cached value of log (2 gamma).
00096         interval log2gamma;
00097 
00098         /// An auxiliary function for the computation of the root of degree
00099         /// gamma of the absolute value of a number.
00100         /// Returns the upper or lower bound for the result.
00101         numType gammaRoot (numType x, bool upperBound) const;
00102 
00103 }; /* mapUnimodalIntv */
00104 
00105 // --------------------------------------------------
00106 
00107 template <class numType>
00108 inline mapUnimodalIntv<numType>::mapUnimodalIntv (const numType &_gamma1,
00109         const numType &_gamma2):
00110         gamma (interval (_gamma1, _gamma1) / interval (_gamma2, _gamma2)),
00111         log2gamma (log (gamma * 2))
00112 {
00113         return;
00114 } /* mapUnimodalIntv::mapUnimodalIntv */
00115 
00116 template <class numType>
00117 inline std::string mapUnimodalIntv<numType>::name () const
00118 {
00119         std::ostringstream nameStr;
00120         nameStr << "unimodal-intv-" << gamma. leftBound ();
00121         return nameStr. str ();
00122 } /* mapUnimodalIntv::name */
00123 
00124 template <class numType>
00125 inline int mapUnimodalIntv<numType>::countCritical () const
00126 {
00127         return 1;
00128 } /* mapUnimodalIntv::countCritical */
00129 
00130 template <class numType>
00131 inline numType mapUnimodalIntv<numType>::criticalPoint (int n) const
00132 {
00133         return 0.0;
00134 } /* mapUnimodalIntv::criticalPoint */
00135 
00136 template <class numType>
00137 inline numType mapUnimodalIntv<numType>::leftBound () const
00138 {
00139         return -1;
00140 } /* mapUnimodalIntv::leftBound */
00141 
00142 template <class numType>
00143 inline numType mapUnimodalIntv<numType>::rightBound () const
00144 {
00145         return 1;
00146 } /* mapUnimodalIntv::rightBound */
00147 
00148 template <class numType>
00149 inline numType mapUnimodalIntv<numType>::gammaRoot (numType x,
00150         bool upperBound) const
00151 {
00152         // return 0 if x is zero
00153         const numType zero (0);
00154         if (x == zero)
00155                 return zero;
00156 
00157         // make x positive if it is negative
00158         if (x < zero)
00159                 x = -x;
00160 
00161         // return the root of degree gamma of x
00162         if (upperBound)
00163                 return exp (log (interval (x, x)) / gamma). rightBound ();
00164         else
00165                 return exp (log (interval (x, x)) / gamma). leftBound ();
00166 
00167 } /* mapUnimodalIntv::gammaRoot */
00168 
00169 template <class numType>
00170 inline void mapUnimodalIntv<numType>::image (const numType &x1,
00171         const numType &x2, numType &y1, numType &y2) const
00172 {
00173         if (x2 < x1)
00174                 throw "Image computation: Wrong interval for 'x'.";
00175         if ((x1 <= 0) && (0 <= x2))
00176         {
00177                 const numType x = (-x1 < x2) ? x2 : -x1;
00178                 y1 = -(this -> paramMax);
00179                 y2 = (power (interval (x, x), gamma) * 2 -
00180                         this -> paramMin). rightBound ();
00181         }
00182         else
00183         {
00184                 interval img = power ((x2 < 0) ? interval (-x2, -x1) :
00185                         interval (x1, x2), gamma) * 2 -
00186                         interval (this -> paramMin, this -> paramMax);
00187                 y1 = img. leftBound ();
00188                 y2 = img. rightBound ();
00189         }
00190         capd::rounding::DoubleRounding::roundNearest ();
00191         return;
00192 } /* mapUnimodalIntv::image */
00193 
00194 template <class numType>
00195 inline numType mapUnimodalIntv<numType>::minLogDerivative (const numType &x1,
00196         const numType &x2, const numType &y1, const numType &y2) const
00197 {
00198         // make sure the input data is correct
00199         if (x2 < x1)
00200                 throw "MinLogDerivative: Wrong interval for 'x'.";
00201         if (y2 < y1)
00202                 throw "MinLogDerivative: Wrong interval for 'y'.";
00203         if (((x1 < 0) || (x1 == 0)) && ((x2 == 0) || (0 < x2)))
00204                 throw "MinLogDerivative: The interval contains zero.";
00205 
00206         // compute the positive preimage of the interval [y1,y2],
00207         // intersect the computed preimage with the interval [x1,x2],
00208         // and return the log of the derivative at the endpoint at which
00209         // the derivative is minimal (right for gamma < 1, left otherwise)
00210         numType result;
00211         if (gamma. rightBound () < 1)
00212         {
00213                 const numType sum = ((interval (y2, y2) +
00214                         this -> paramMax) / 2). rightBound ();
00215                 const numType preImg = (0 < sum) ?
00216                         gammaRoot (sum, true) : 0;
00217                 const numType x = (0 < x1) ? x2 : -x1;
00218                 const numType endPoint = (x < preImg) ? x : preImg;
00219                 result = (log2gamma + (gamma - 1) *
00220                         log (interval (endPoint, endPoint))). leftBound ();
00221         }
00222         else
00223         {
00224                 const numType sum = ((interval (y1, y1) +
00225                         this -> paramMin) / 2). leftBound ();
00226                 if ((sum < 0) || (sum == 0))
00227                         throw "MinLogDerivative: Zero in the preimage.";
00228                 const numType preImg = (0 < sum) ?
00229                         gammaRoot (sum, false) : 0;
00230                 const numType x = (0 < x1) ? x1 : -x2;
00231                 const numType endPoint = (x < preImg) ? preImg : x;
00232                 result = (log2gamma + (gamma - 1) *
00233                         log (interval (endPoint, endPoint))). leftBound ();
00234         }
00235         capd::rounding::DoubleRounding::roundNearest ();
00236         return result;
00237 } /* mapUnimodalIntv::minLogDerivative */
00238 
00239 
00240 } // namespace unifexp
00241 
00242 #endif // _mapunimi_h_
00243 
00244 /// @}
00245 

Generated on Wed Nov 21 11:08:41 2007 for The Uniform Expansion Software by  doxygen 1.5.3