Main Page   Namespace List   Class Hierarchy   Alphabetical List   Compound List   File List   Namespace Members   Compound Members   File Members   Related Pages  

odesystem.h

00001 /***************************************************************************
00002                           odesystem.h  -  description
00003                              -------------------
00004     begin                : Tue Jul 11 2000
00005     copyright            : (C) 2000 by Michael Peeters
00006     email                : Michael.Peeters@vub.ac.be
00007  ***************************************************************************/
00008 
00009 /***************************************************************************
00010  *                                                                         *
00011  *   This program is free software; you can redistribute it and/or modify  *
00012  *   it under the terms of the GNU General Public License as published by  *
00013  *   the Free Software Foundation; either version 2 of the License, or     *
00014  *   (at your option) any later version.                                   *
00015  *                                                                         *
00016  ***************************************************************************/
00017 
00018 #ifndef ODESYSTEM_H
00019 #define ODESYSTEM_H
00020 
00021 #include "rootscan.h"
00022 
00023 #include "timeframe.h"
00024 #include "random.h"
00025 // #include "jacobian.h"
00026 // #include "debugmacro.h"
00027 #include "ticktock.h"
00028 #include "integrator.h"
00029 #include "cowner.h"
00030 #include "cycler.h"
00031 #include "numerictraits.h"
00032 #include "numerictypes.h"
00033 
00034 namespace MODEL {
00035 
00048   template<integer dims, typename nelem=number, class NT = NumericTraits<nelem,dims> >
00049   class ODESystem : public TickTock
00050   {
00051   public:
00052         typedef typename NT::number                                     numT;
00053         typedef typename NT::vect                                       vect;
00054         typedef typename NT::vf                                         vf;
00055         typedef typename NT::matrix                                     matrix;
00056         
00057         
00062         ODESystem(      TimeFrame& T, vf& detpart, vf& stochpart, vect& initial,
00063                                 Integration it=Euler, time res=1E-3);
00064 
00067         ODESystem(      TimeFrame& T, vect& initial, Integrator<dims,nelem,NT>* i, time res=1e-3);
00068 
00070         ODESystem(      TimeFrame& T, vf& detpart, vect& initial,
00071                                 Integration it, time res);
00072 
00073         virtual ~ODESystem();
00074 
00077         virtual void tick()
00078         {
00079           step(); // Take care of the TimeFrame aspect (the modulators...)
00080           integ->step(current,dt); // Call the integrator
00081         }
00082 
00084         const vect& get_current(void)
00085         {
00086           return current;
00087         }
00088 
00089         void set_current(const vect& hereandnow)
00090         {
00091           current = hereandnow;  
00092         }
00093 
00098         const vect&     operator()(void){return get_current();}
00099 
00101         void reset(void) 
00102         {
00103           reset(init);
00104         }
00105         
00106         
00108         void    reset(const vect& fromhere)
00109         {  
00110           current=fromhere; 
00111           //Previously on MODEL:
00112           //relax_independent();
00113           
00114           // But a much beter way to go is:
00115           RootScan<dims,nelem,NT> find(*deter);
00116           find.add_start(fromhere);
00117           ScanList<dims,nelem,NT> roots=find.scan(0.,1.,1);
00118           ScanList<dims,nelem,NT> goodone=roots.select(RootScan<dims>::stable);
00119           
00120           if (goodone.get_data()[0.].size()>0)
00121                 {         
00122                   current=goodone.get_data()[0.].front()[0];
00123                   cerr << "ODESystem: started at " << current << endl;
00124                 }
00125           else
00126                 {
00127                   current=fromhere;
00128                   cerr << "ODEsystem - No stable starting point found @ " <<
00129                         current << endl;
00130                 }
00131         }
00132 
00136         time relax(number goal=DEFAULT_RELAX);
00139         time relax_independent(number goal=DEFAULT_RELAX);
00140         
00143   private:
00144         // Variables
00145         static const number DEFAULT_RELAX=1E-7;
00146 
00148         vect            init;
00149 
00151         vect current;
00152 
00154         Integrator<dims,nelem,NT>*      integ;
00155 
00157         vf* deter;
00159         vf* stoch;
00160   };
00161 
00162   // Member Functions
00163   //--------------------------------------------------------------------------------
00164 
00165   template<integer dims, typename nelem, class NT >
00166   ODESystem<dims,nelem,NT>::ODESystem(  TimeFrame& T, vf& detpart, vf& stochpart, vect& initial,
00167                                                                                 Integration it, time res) :
00168         TickTock(T,res), init(initial), integ(NULL), deter(&detpart), stoch(&stochpart)
00169   {
00170         // dt might have changed (because it has to fit integrally inside the top
00171         // TimeFrame)
00172 
00173         // This is not needed: dt is a member !!!!
00174         // dt=get_dt();
00175         // relax to stable state to start if possible
00176 
00177         current = init;
00178         reset();
00179         
00180         // Initialize integrator
00181         switch(it)
00182           {
00183           case Euler:   integ = new IEuler<dims,nelem,NT>(detpart,stochpart);
00184                 break;
00185           case RungeKutta: integ = new IRungeKutta<dims,nelem,NT>(detpart,stochpart);
00186                 break;
00187           case Milshtein: integ = new IMilshtein<dims,nelem,NT>(detpart,stochpart);
00188                 break;
00189 
00190           case Heun: integ = new IHeun<dims,nelem,NT>(detpart,stochpart);
00191                 break;
00192           
00193           default:      throw std::logic_error("ODESystem::Integration method not implemented");
00194                 break;
00195           }
00196   }
00197   //-------------------------------------------------------------------------------
00198 
00201   template<integer dims, typename nelem, class NT >
00202         ODESystem<dims,nelem,NT>::ODESystem(    TimeFrame& T, vect& initial, Integrator<dims,nelem,NT>* i, time res=1e-3):
00203         TickTock(T,res), init(initial), integ(NULL)
00204   {
00205         // dt might have changed (because it has to fit integrally inside the top
00206         // TimeFrame)
00207 
00208         // This is not needed: dt is a member !!!!
00209         // dt=get_dt();
00210         // relax to stable state to start if possible
00211 
00212         integ = i;
00213 
00214         // get functions !!!
00215         deter = i->get_deterministic() ; stoch= i->get_stochastic();
00216 
00217         current = init;
00218         reset();        
00219   }
00220 
00221   //--------------------------------------------------------------------------------
00222 
00225   template<integer dims, typename nelem, class NT >
00226   ODESystem<dims,nelem,NT>::ODESystem(  TimeFrame& T, vf& detpart, vect& initial,
00227                                                                                 Integration it, time res) :
00228         TickTock(T,res), init(initial), integ(NULL), deter(&detpart)
00229   {
00230         // dt might have changed (because it has to fit integrally inside the top
00231         // TimeFrame)
00232 
00233         // This is not needed: dt is a member !!!!
00234         // dt=get_dt();
00235         // relax to stable state to start if possible
00236 
00237         current = init;
00238         reset();
00239         
00240         // Initialize integrator
00241         switch(it)
00242           {
00243           case Euler:   integ = new IEuler<dims,nelem,NT>(detpart,detpart);
00244                 break;
00245           case RungeKutta: integ = new IRungeKutta<dims,nelem,NT>(detpart,detpart);
00246                 break;
00247           
00248           default:      throw std::logic_error("ODESystem::Integration method not implemented");
00249                 break;
00250           }
00251   }
00252 
00253   //--------------------------------------------------------------------------------
00254 
00255   template<integer dims, typename nelem, class NT >
00256   ODESystem<dims,nelem,NT>::~ODESystem()
00257   {
00258         delete integ;
00259   }
00260 
00261   //--------------------------------------------------------------------------------
00262 
00263   template<integer dims, typename nelem, class NT >
00264   time 
00265   ODESystem<dims,nelem,NT>::relax_independent(number goal)
00266   {
00267         // This will become a RungeKutta
00268         IRungeKutta<dims,nelem,NT> massage(*deter,*stoch); 
00269         ODESystem& system(*this); // alias for readability
00270         number change(0.);
00271         time chrono=0.;
00272         vect lastvalue(system());
00273 
00274         do 
00275           {
00276                 // One step
00277                 chrono+=dt;
00278                 vect newvalue=massage.step(current,dt);
00279                                   
00280                 // Have we converged?
00281                 change=0.;
00282                 for (counter i=0;i<dims;i++)
00283                   {
00284                         number converge=abs((newvalue[i]-lastvalue[i])/lastvalue[i]);
00285                         if(converge>change) change=converge; // find maximal change
00286                   }
00287                 lastvalue=newvalue;
00288                 // cout << change << "\t" << newvalue << endl;
00289         
00290           } while (change>goal);
00291 
00292         return chrono; // return time spent
00293   }
00294 
00295   //--------------------------------------------------------------------------------
00296   template<integer dims, typename nelem, class NT >
00297   time 
00298   ODESystem<dims,nelem,NT>::relax(number goal)
00299   {
00300         ODESystem& system(*this); // alias for readability
00301         number change(0.);
00302         TimeFrame& t=get_timeframe();
00303         time init=t;
00304         vect converge(0.),lastvalue(system());
00305 
00306         do 
00307           {
00308                 // One step
00309                 ++t;
00310                 vect newvalue=system();
00311                   
00312                 // Have we converged?
00313                 change=0.;
00314                 for (counter i=0;i<dims;i++)
00315                   {
00316                         number converge=abs((newvalue[i]-lastvalue[i])/lastvalue[i]);
00317                         if(converge>change) change=converge; // find maximal change
00318                   }
00319                 lastvalue=newvalue;
00320 
00321                 
00322                 //              DEBUGSTREAM << change <<endl;
00323           } while (change>goal);
00324 
00325         return t-init; // return time spent
00326   }
00327 
00328 } // end namespace
00329 #endif
00330 
00331 /*********************************************************************
00332 $Id: odesystem.h,v 1.2.2.2 2002/06/19 13:17:32 mpeeters Exp $
00333 **********************************************************************
00334 
00335 $Log: odesystem.h,v $
00336 Revision 1.2.2.2  2002/06/19 13:17:32  mpeeters
00337 Possibility for custom integrator.
00338 
00339 Revision 1.2.2.1  2001/10/15 15:08:24  mpeeters
00340 Added staring from a different point to make restarting easier.
00341 
00342 Revision 1.2  2001/07/26 11:53:28  mpeeters
00343 Added constructor for non-stochastic case.
00344 Made sure the initial stat point is found correctly.
00345 Removed unneeded get_dt().
00346 
00347 Revision 1.1  2001/05/22 10:54:55  mpeeters
00348 Moved sources and headers for libModel to model/ subdirectory, in an attempt to rationalize the source tree. This should make things "netter".
00349 
00350 Revision 1.5  2001/05/21 11:53:16  mpeeters
00351 Removed Makefile.in, which is automatically generated anyway.
00352 
00353 Revision 1.4  2000/09/29 08:55:07  mpeeters
00354 Inserted check for stable state. Removed bug where system would choose wrong starting point.
00355 
00356 Revision 1.3  2000/09/22 08:56:04  mpeeters
00357 Bug fix: Rootscan was called for a 3D system always. Changed template
00358 parameter to RootScan<dims>, as it should be.
00359 
00360 Revision 1.2  2000/09/15 10:26:31  mpeeters
00361 Cleaned out header and added CVS tails to files, removed superfuous
00362 @author comments, inserted dates
00363 
00364 *********************************************************************/

To get the sources or tarballs, please go to SourceForge or you can use the CVS repository.

More Info? Michael Peeters. Also, check our research website: www.alna.vub.ac.be

Last update: June 2002.


Looking for Open Source? Check out SourceForge Logo !