mapcubi.h

Go to the documentation of this file.
00001 /// @addtogroup unifexp
00002 /// @{
00003 
00004 /////////////////////////////////////////////////////////////////////////////
00005 ///
00006 /// @file mapcubi.h
00007 ///
00008 /// This file contains the definition of the cubic map -x^3 + 3x - a
00009 /// on [-4,4] 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 24, 2007. Last revision: August 25, 2007.
00034 
00035 #ifndef _mapcubi_h_
00036 #define _mapcubi_h_
00037 
00038 #include <string>
00039 #include <iostream>
00040 #include "maptype.h"
00041 
00042 // the "interval" package from the CAPD library
00043 #include "interval/doubleInterval.h"
00044 
00045 
00046 namespace unifexp {
00047 
00048 // --------------------------------------------------
00049 // ----------------- the cubic map ------------------
00050 // --------------------------------------------------
00051 
00052 /// This class defines the cubic map -x^3 + 3x - a on [-4,4]
00053 /// with the usage of interval arithmetic.
00054 /// It is suitable for rigorous computations.
00055 /// Good values of the parameter are between 0 and 2,
00056 /// more or less (the optimal value is 2, I think).
00057 template <class numType>
00058 class mapCubicIntv: public mapType<numType>
00059 {
00060 public:
00061         /// The constructor.
00062         mapCubicIntv ();
00063 
00064         /// Returns the name of the object.
00065         std::string name () const;
00066 
00067         /// Returns the number of critical points.
00068         int countCritical () const;
00069 
00070         /// Returns the subsequent critical points.
00071         numType criticalPoint (int n) const;
00072 
00073         /// Returns the left bound of the domain of the map.
00074         numType leftBound () const;
00075 
00076         /// Returns the right bound of the domain of the map.
00077         numType rightBound () const;
00078 
00079         /// Computes an enclosure of the image of the given interval.
00080         void image (const numType &x1, const numType &x2,
00081                 numType &y1, numType &y2) const;
00082 
00083         /// Computes the minimal log of the derivative over those points
00084         /// in the interval [x1,x2] whose images may fall into [y1,y2]
00085         numType minLogDerivative (const numType &x1, const numType &x2,
00086                 const numType &y1, const numType &y2) const;
00087 
00088 }; /* mapCubicIntv */
00089 
00090 // --------------------------------------------------
00091 
00092 template <class numType>
00093 inline mapCubicIntv<numType>::mapCubicIntv ()
00094 {
00095         return;
00096 } /* mapCubicIntv::mapCubicIntv */
00097 
00098 template <class numType>
00099 inline std::string mapCubicIntv<numType>::name () const
00100 {
00101         return std::string ("cubic-intv");
00102 } /* mapCubicIntv::name */
00103 
00104 template <class numType>
00105 inline int mapCubicIntv<numType>::countCritical () const
00106 {
00107         return 2;
00108 } /* mapCubicIntv::countCritical */
00109 
00110 template <class numType>
00111 inline numType mapCubicIntv<numType>::criticalPoint (int n) const
00112 {
00113         return n ? 1 : -1;
00114 } /* mapCubicIntv::criticalPoint */
00115 
00116 template <class numType>
00117 inline numType mapCubicIntv<numType>::leftBound () const
00118 {
00119         return -4;
00120 } /* mapCubicIntv::leftBound */
00121 
00122 template <class numType>
00123 inline numType mapCubicIntv<numType>::rightBound () const
00124 {
00125         return 4;
00126 } /* mapCubicIntv::rightBound */
00127 
00128 template <class numType>
00129 inline void mapCubicIntv<numType>::image (const numType &x1,
00130         const numType &x2, numType &y1, numType &y2) const
00131 {
00132         if (x2 < x1)
00133                 throw "Image computation: Wrong interval for 'x'.";
00134         interval x1i (x1, x1);
00135         interval x2i (x2, x2);
00136         if ((x2 <= -1) || (1 <= x1))
00137         {
00138                 y1 = (-(x2i * x2i - 3) * x2i - this -> paramMax).
00139                         leftBound ();
00140                 y2 = (-(x1i * x1i - 3) * x1i - this -> paramMin).
00141                         rightBound ();
00142         }
00143         else if ((-1 <= x1) && (x2 <= 1))
00144         {
00145                 y1 = (-(x1i * x1i - 3) * x1i - this -> paramMax).
00146                         leftBound ();
00147                 y2 = (-(x2i * x2i - 3) * x2i - this -> paramMin).
00148                         rightBound ();
00149         }
00150         else if (x2 <= 1)
00151         {
00152                 const interval y3 = -(x1i * x1i - 3) * x1i;
00153                 const interval y4 = -(x2i * x2i - 3) * x2i;
00154                 y1 = (interval (-2, -2) - this -> paramMax).
00155                         leftBound ();
00156                 y2 = (((y3. rightBound () < y4. rightBound ()) ? y4 : y3) -
00157                         this -> paramMin). rightBound ();
00158         }
00159         else if (-1 <= x1)
00160         {
00161                 const interval y3 = -(x1i * x1i - 3) * x1i;
00162                 const interval y4 = -(x2i * x2i - 3) * x2i;
00163                 y1 = (((y3. leftBound () < y4. leftBound ()) ? y3 : y4) -
00164                         this -> paramMax). leftBound ();
00165                 y2 = (interval (4, 4) - this -> paramMin). rightBound ();
00166         }
00167         else
00168                 throw "Cubic map: Too large domain interval.";
00169         capd::rounding::DoubleRounding::roundNearest ();
00170         return;
00171 } /* mapCubicIntv::image */
00172 
00173 template <class numType>
00174 inline numType mapCubicIntv<numType>::minLogDerivative (const numType &x1,
00175         const numType &x2, const numType &y1, const numType &y2) const
00176 {
00177         // make sure the input data is correct
00178         if (x2 < x1)
00179                 throw "MinLogDerivative: Wrong interval for 'x'.";
00180         if (y2 < y1)
00181                 throw "MinLogDerivative: Wrong interval for 'y'.";
00182 
00183         // NOTE: Since it is not so easy to compute the preimage of the
00184         // interval [y1,y2] with respect to the cubic polynomial map,
00185         // the lower bound for the log of the derivative is computed
00186         // solely based on the interval [x1,x2]
00187         interval x1i (x1, x1);
00188         interval x2i (x2, x2);
00189         numType result;
00190         if (x2 < -1)
00191         {
00192                 result = log (3 * x2i * x2i + 3). leftBound ();
00193         }
00194         else if (1 < x1)
00195         {
00196                 result = log (3 * x1i * x1i + 3). leftBound ();
00197         }
00198         else if ((-1 < x1) && (x2 < 1))
00199         {
00200                 if (0 < x1)
00201                         result = log (-3 * x2i * x2i + 3). leftBound ();
00202                 else if (x2 < 0)
00203                         result = log (-3 * x1i * x1i + 3). leftBound ();
00204                 else
00205                 {
00206                         result = log (-3 * ((x1 < x2) ? (x2i * x2i) :
00207                                 (x1i * x1i)) + 3). leftBound ();
00208                 }
00209         }
00210         else
00211         {
00212                 throw "Log of interval containing zero.";
00213                 result = 0;
00214         }
00215         capd::rounding::DoubleRounding::roundNearest ();
00216         return result;
00217 } /* mapCubicIntv::minLogDerivative */
00218 
00219 
00220 } // namespace unifexp
00221 
00222 #endif // _mapcubi_h_
00223 
00224 /// @}
00225 

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