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 *********************************************************************/
More Info? Michael Peeters. Also, check our research website: www.alna.vub.ac.be
Last update: June 2002.
Looking for Open Source? Check out
!