DifferentialStepper.hpp

Go to the documentation of this file.
00001 //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
00002 //
00003 //        This file is part of E-Cell Simulation Environment package
00004 //
00005 //                Copyright (C) 1996-2002 Keio University
00006 //
00007 //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
00008 //
00009 //
00010 // E-Cell is free software; you can redistribute it and/or
00011 // modify it under the terms of the GNU General Public
00012 // License as published by the Free Software Foundation; either
00013 // version 2 of the License, or (at your option) any later version.
00014 // 
00015 // E-Cell is distributed in the hope that it will be useful,
00016 // but WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
00018 // See the GNU General Public License for more details.
00019 // 
00020 // You should have received a copy of the GNU General Public
00021 // License along with E-Cell -- see the file COPYING.
00022 // If not, write to the Free Software Foundation, Inc.,
00023 // 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
00024 // 
00025 //END_HEADER
00026 //
00027 // written by Koichi Takahashi <shafi@e-cell.org>,
00028 // E-Cell Project.
00029 //
00030 
00031 #ifndef __DIFFERENTIALSTEPPER_HPP
00032 #define __DIFFERENTIALSTEPPER_HPP
00033 
00034 #include "libecs.hpp"
00035 #include "Stepper.hpp"
00036 
00037 #include "boost/multi_array.hpp"
00038 
00039 namespace libecs
00040 {
00041 
00042   /** @addtogroup stepper
00043    *@{
00044    */
00045 
00046   /** @file */
00047 
00048   /**
00049      DIFFERENTIAL EQUATION SOLVER
00050 
00051 
00052   */
00053 
00054   //  DECLARE_VECTOR( RealVector, RealMatrix );
00055   
00056   typedef boost::multi_array<Real, 2> RealMatrix_;
00057   DECLARE_TYPE( RealMatrix_, RealMatrix );
00058 
00059 
00060   DECLARE_CLASS( DifferentialStepper );
00061 
00062   LIBECS_DM_CLASS( DifferentialStepper, Stepper )
00063   {
00064 
00065   public:
00066 
00067     LIBECS_DM_OBJECT_ABSTRACT( DifferentialStepper )
00068       {
00069         INHERIT_PROPERTIES( Stepper );
00070 
00071         // FIXME: load/save ??
00072         PROPERTYSLOT( Real, StepInterval,
00073                       &DifferentialStepper::initializeStepInterval,
00074                       &DifferentialStepper::getStepInterval );
00075         
00076         PROPERTYSLOT_GET_NO_LOAD_SAVE( Real, NextStepInterval );
00077         PROPERTYSLOT_SET_GET_NO_LOAD_SAVE( Real,  TolerableStepInterval );
00078         PROPERTYSLOT_GET_NO_LOAD_SAVE( Integer,  Stage );
00079         PROPERTYSLOT_GET_NO_LOAD_SAVE( Integer,  Order );
00080       }
00081 
00082     class Interpolant
00083       :
00084       public libecs::Interpolant
00085     {
00086 
00087     public:
00088 
00089       Interpolant( DifferentialStepperRef aStepper, 
00090                      VariablePtr const aVariablePtr )
00091         :
00092         libecs::Interpolant( aVariablePtr ),
00093         theStepper( aStepper ),
00094         theIndex( theStepper.getVariableIndex( aVariablePtr ) )
00095       {
00096         ; // do nothing
00097       }
00098       
00099 
00100       /*
00101         The getDifference() below is an optimized version of
00102         the original implementation based on the following two functions.
00103         (2004/10/19)
00104 
00105       const Real interpolate( const RealMatrixCref aTaylorSeries,
00106                               const Real anInterval,
00107                               const Real aStepInterval )
00108       {
00109         const Real theta( anInterval / aStepInterval );
00110 
00111         Real aDifference( 0.0 );
00112         Real aFactorialInv( 1.0 );
00113 
00114         for ( RealMatrix::size_type s( 0 ); s < aTaylorSeries.size(); ++s )
00115           {
00116             //      aFactorialInv /= s + 1;
00117             aDifference += aTaylorSeries[ s ][ theIndex ] * aFactorialInv;
00118             aFactorialInv *= theta;
00119           }
00120 
00121         return aDifference * anInterval;
00122       }
00123 
00124       virtual const Real getDifference( RealParam aTime, 
00125                                         RealParam anInterval )
00126       {
00127 
00128         if ( !theStepper.theStateFlag )
00129           {
00130             return 0.0;
00131           }
00132 
00133         const RealMatrixCref aTaylorSeries( theStepper.getTaylorSeries() );
00134         const Real aTimeInterval( aTime - theStepper.getCurrentTime() );
00135         const Real aStepInterval( theStepper.getTolerableStepInterval() );
00136         
00137 
00138         const Real i1( interpolate( aTaylorSeries, 
00139                                     aTimeInterval, 
00140                                     aStepInterval ) );
00141         const Real i2( interpolate( aTaylorSeries, 
00142                                     aTimeInterval - anInterval, 
00143                                     aStepInterval ) );
00144         return ( i1 - i2 );
00145 
00146         }
00147       */
00148 
00149       virtual const Real getDifference( RealParam aTime, 
00150                                         RealParam anInterval ) const
00151       {
00152         if ( !theStepper.theStateFlag )
00153           {
00154             return 0.0;
00155           }
00156 
00157         const Real aTimeInterval1( aTime - theStepper.getCurrentTime() );
00158         const Real aTimeInterval2( aTimeInterval1 - anInterval );
00159 
00160         RealMatrixCref aTaylorSeries( theStepper.getTaylorSeries() );
00161         RealCptr aTaylorCoefficientPtr( aTaylorSeries.origin() + theIndex );
00162 
00163         // calculate first order.
00164         // here it assumes that always aTaylorSeries.size() >= 1
00165 
00166         // *aTaylorCoefficientPtr := aTaylorSeries[ 0 ][ theIndex ]
00167         Real aValue1( *aTaylorCoefficientPtr * aTimeInterval1 );
00168         Real aValue2( *aTaylorCoefficientPtr * aTimeInterval2 );
00169 
00170 
00171         // check if second and higher order calculations are necessary.
00172         //      const RealMatrix::size_type aTaylorSize( aTaylorSeries.size() );
00173         const RealMatrix::size_type aTaylorSize( theStepper.getOrder() );
00174         if( aTaylorSize >= 2)
00175           {
00176             const Real 
00177               aStepIntervalInv( 1.0 / theStepper.getTolerableStepInterval() );
00178             
00179             const RealMatrix::size_type aStride( aTaylorSeries.strides()[0] );
00180 
00181             Real aFactorialInv1( aTimeInterval1 );
00182             Real aFactorialInv2( aTimeInterval2 );
00183 
00184             RealMatrix::size_type s( aTaylorSize - 1 );
00185 
00186             const Real theta1( aTimeInterval1 * aStepIntervalInv );
00187             const Real theta2( aTimeInterval2 * aStepIntervalInv );
00188 
00189             do 
00190               {
00191                 // main calculation for the 2+ order
00192                 
00193                 // aTaylorSeries[ s ][ theIndex ]
00194                 aTaylorCoefficientPtr += aStride;
00195                 const Real aTaylorCoefficient( *aTaylorCoefficientPtr );
00196                 
00197                 aFactorialInv1 *= theta1;
00198                 aFactorialInv2 *= theta2;
00199                 
00200                 aValue1 += aTaylorCoefficient * aFactorialInv1;
00201                 aValue2 += aTaylorCoefficient * aFactorialInv2;
00202                 
00203                 // LIBECS_PREFETCH( aTaylorCoefficientPtr + aStride, 0, 1 );
00204                 --s;
00205               } while( s != 0 );
00206           }
00207 
00208         return aValue1 - aValue2;
00209       }
00210       
00211       virtual const Real getVelocity( RealParam aTime ) const
00212       {
00213         if ( !theStepper.theStateFlag )
00214           {
00215             return 0.0;
00216           }
00217 
00218         const Real aTimeInterval( aTime - theStepper.getCurrentTime() );
00219 
00220         RealMatrixCref aTaylorSeries( theStepper.getTaylorSeries() );
00221         RealCptr aTaylorCoefficientPtr( aTaylorSeries.origin() + theIndex );
00222 
00223         // calculate first order.
00224         // here it assumes that always aTaylorSeries.size() >= 1
00225 
00226         // *aTaylorCoefficientPtr := aTaylorSeries[ 0 ][ theIndex ]
00227         Real aValue( *aTaylorCoefficientPtr );
00228 
00229         // check if second and higher order calculations are necessary.
00230         //      const RealMatrix::size_type aTaylorSize( aTaylorSeries.size() );
00231 
00232         const RealMatrix::size_type aTaylorSize( theStepper.getStage() );
00233         if( aTaylorSize >= 2 && aTimeInterval != 0.0 )
00234           {
00235             const RealMatrix::size_type aStride( aTaylorSeries.strides()[0] );
00236 
00237             Real aFactorialInv( 1.0 );
00238 
00239             RealMatrix::size_type s( 1 );
00240 
00241             const Real theta( aTimeInterval 
00242                               / theStepper.getTolerableStepInterval() );
00243 
00244             do 
00245               {
00246                 // main calculation for the 2+ order
00247                 ++s;
00248                 
00249                 aTaylorCoefficientPtr += aStride;
00250                 const Real aTaylorCoefficient( *aTaylorCoefficientPtr );
00251                 
00252                 aFactorialInv *= theta * s;
00253                 
00254                 aValue += aTaylorCoefficient * aFactorialInv;
00255                 
00256                 // LIBECS_PREFETCH( aTaylorCoefficientPtr + aStride, 0, 1 );
00257               } while( s != aTaylorSize );
00258           }
00259 
00260         return aValue;
00261       }
00262       
00263     protected:
00264 
00265       DifferentialStepperRef    theStepper;
00266       VariableVector::size_type theIndex;
00267 
00268     };
00269 
00270   public:
00271 
00272     ECELL_API DifferentialStepper();
00273 
00274     ECELL_API virtual ~DifferentialStepper();
00275 
00276     SET_METHOD( Real, NextStepInterval )
00277     {
00278       theNextStepInterval = value;
00279     }
00280 
00281     GET_METHOD( Real, NextStepInterval )
00282     {
00283       return theNextStepInterval;
00284     }
00285 
00286     SET_METHOD( Real, TolerableStepInterval )
00287     {
00288       theTolerableStepInterval = value;
00289     }
00290 
00291     GET_METHOD( Real, TolerableStepInterval )
00292     {
00293       return theTolerableStepInterval;
00294     }
00295 
00296     void initializeStepInterval( RealParam aStepInterval )
00297     {
00298       setStepInterval( aStepInterval );
00299       setTolerableStepInterval( aStepInterval );
00300       setNextStepInterval( aStepInterval );
00301     }
00302 
00303     ECELL_API void resetAll();
00304 
00305     ECELL_API void interIntegrate();
00306 
00307     void initializeVariableReferenceList();
00308 
00309     ECELL_API void setVariableVelocity( boost::detail::multi_array::sub_array<Real, 1> aVelocityBuffer );
00310  
00311     ECELL_API virtual void initialize();
00312 
00313     ECELL_API virtual void reset();
00314 
00315     ECELL_API virtual void interrupt( TimeParam aTime );
00316 
00317     virtual InterpolantPtr createInterpolant( VariablePtr aVariable )
00318     {
00319       return new DifferentialStepper::Interpolant( *this, aVariable );
00320     }
00321 
00322     virtual GET_METHOD( Integer, Stage )
00323     { 
00324       return 1; 
00325     }
00326 
00327     virtual GET_METHOD( Integer, Order )
00328     { 
00329       return getStage(); 
00330     }
00331 
00332     RealMatrixCref getTaylorSeries() const
00333     {
00334       return theTaylorSeries;
00335     }
00336 
00337   protected:
00338 
00339     const bool isExternalErrorTolerable() const;
00340 
00341     RealMatrix theTaylorSeries;
00342 
00343     std::vector< std::vector<Integer> > theVariableReferenceListVector;
00344 
00345     bool theStateFlag;
00346 
00347   private:
00348 
00349     Real theNextStepInterval;
00350     Real theTolerableStepInterval;
00351   };
00352 
00353 
00354   /**
00355      ADAPTIVE STEPSIZE DIFFERENTIAL EQUATION SOLVER
00356 
00357 
00358   */
00359 
00360   DECLARE_CLASS( AdaptiveDifferentialStepper );
00361 
00362   LIBECS_DM_CLASS( AdaptiveDifferentialStepper, DifferentialStepper )
00363   {
00364 
00365   public:
00366 
00367     LIBECS_DM_OBJECT_ABSTRACT( AdaptiveDifferentialStepper )
00368       {
00369         INHERIT_PROPERTIES( DifferentialStepper );
00370 
00371         PROPERTYSLOT_SET_GET( Real, Tolerance );
00372         PROPERTYSLOT_SET_GET( Real, AbsoluteToleranceFactor );
00373         PROPERTYSLOT_SET_GET( Real, StateToleranceFactor );
00374         PROPERTYSLOT_SET_GET( Real, DerivativeToleranceFactor );
00375 
00376         PROPERTYSLOT( Integer, IsEpsilonChecked,
00377                       &AdaptiveDifferentialStepper::setEpsilonChecked,
00378                       &AdaptiveDifferentialStepper::isEpsilonChecked );
00379         PROPERTYSLOT_SET_GET( Real, AbsoluteEpsilon );
00380         PROPERTYSLOT_SET_GET( Real, RelativeEpsilon );
00381 
00382         PROPERTYSLOT_GET_NO_LOAD_SAVE( Real, MaxErrorRatio );
00383       }
00384 
00385   public:
00386 
00387     ECELL_API AdaptiveDifferentialStepper();
00388 
00389     ECELL_API virtual ~AdaptiveDifferentialStepper();
00390 
00391     /**
00392        Adaptive stepsize control.
00393 
00394        These methods are for handling the standerd error control objects.
00395     */
00396 
00397     SET_METHOD( Real, Tolerance )
00398     {
00399       theTolerance = value;
00400     }
00401 
00402     GET_METHOD( Real, Tolerance )
00403     {
00404       return theTolerance;
00405     }
00406 
00407     SET_METHOD( Real, AbsoluteToleranceFactor )
00408     {
00409       theAbsoluteToleranceFactor = value;
00410     }
00411 
00412     GET_METHOD( Real, AbsoluteToleranceFactor )
00413     {
00414       return theAbsoluteToleranceFactor;
00415     }
00416 
00417     SET_METHOD( Real, StateToleranceFactor )
00418     {
00419       theStateToleranceFactor = value;
00420     }
00421 
00422     GET_METHOD( Real, StateToleranceFactor )
00423     {
00424       return theStateToleranceFactor;
00425     }
00426 
00427     SET_METHOD( Real, DerivativeToleranceFactor )
00428     {
00429       theDerivativeToleranceFactor = value;
00430     }
00431 
00432     GET_METHOD( Real, DerivativeToleranceFactor )
00433     {
00434       return theDerivativeToleranceFactor;
00435     }
00436 
00437     SET_METHOD( Real, MaxErrorRatio )
00438     {
00439       theMaxErrorRatio = value;
00440     }
00441 
00442     GET_METHOD( Real, MaxErrorRatio )
00443     {
00444       return theMaxErrorRatio;
00445     }
00446 
00447     /**
00448        check difference in one step
00449     */
00450 
00451     SET_METHOD( Integer, EpsilonChecked )
00452     {
00453       if ( value > 0 ) {
00454         theEpsilonChecked = true;
00455       }
00456       else {
00457         theEpsilonChecked = false;
00458       }
00459     }
00460 
00461     const Integer isEpsilonChecked() const
00462     {
00463       return theEpsilonChecked;
00464     }
00465 
00466     SET_METHOD( Real, AbsoluteEpsilon )
00467     {
00468       theAbsoluteEpsilon = value;
00469     }
00470 
00471     GET_METHOD( Real, AbsoluteEpsilon )
00472     {
00473       return theAbsoluteEpsilon;
00474     }
00475 
00476     SET_METHOD( Real, RelativeEpsilon )
00477     {
00478       theRelativeEpsilon = value;
00479     }
00480 
00481     GET_METHOD( Real, RelativeEpsilon )
00482     {
00483       return theRelativeEpsilon;
00484     }
00485 
00486     ECELL_API virtual void initialize();
00487 
00488     ECELL_API virtual void step();
00489 
00490     virtual bool calculate() = 0;
00491 
00492     virtual GET_METHOD( Integer, Stage )
00493     { 
00494       return 2;
00495     }
00496 
00497   private:
00498 
00499     Real safety;
00500     Real theTolerance;
00501     Real theAbsoluteToleranceFactor;
00502     Real theStateToleranceFactor;
00503     Real theDerivativeToleranceFactor;
00504 
00505     bool theEpsilonChecked;
00506     Real theAbsoluteEpsilon;
00507     Real theRelativeEpsilon;
00508 
00509     Real theMaxErrorRatio;
00510   };
00511 
00512 
00513 } // namespace libecs
00514 
00515 #endif /* __DIFFERENTIALSTEPPER_HPP */
00516 
00517 
00518 /*
00519   Do not modify
00520   $Author: sachiboo $
00521   $Revision: 2618 $
00522   $Date: 2006-11-24 13:08:04 +0100 (Fri, 24 Nov 2006) $
00523   $Locker$
00524 */

Generated on Mon Dec 18 07:29:46 2006 for E-CELL C++ libraries (libecs and libemc) 3.1.105 by  doxygen 1.5.1