pBook‎ > ‎Solver:: Integrator:: C++‎ > ‎

### C++11

posted Nov 29, 2014, 5:22 AM by Javad Taghia   [ updated Nov 29, 2014, 5:22 AM ]
https://github.com/michaeloriordan/integrator

 #include "Integrator.hpp" #include "VecDef.hpp" #include #include void fillVec(vec_t& vec, std::function f, double dx) { const int N = vec.size(); for (int i = 0; i < N; ++i) { vec[i] = f(i * dx); } } int main() { double a = 0; double b = 2; int N = 1000; double dx = (b - a) / N; vec_t vec(N, 0.0); auto f = [](double x){return x*x;}; fillVec(vec, f, dx); double sum = integrator::integrate(vec, dx); std::cout << sum << std::endl; return 0; }

 #ifndef __Integrator_ #define __Integrator_ #include #include #include #include #include namespace si_lib { /** The calendar queue class stores a (time,event_name,recurrence_interval) tuple. It is essentially a priority queue whose top priority is always increasing. */ class CalendarQueue { private: /** This class stores an event. An event is represented by a floating-point time and a string identifying what the event is. The event may recur at a regular interval. */ class DiscreteEvent { public: double t; std::string event; double recur_int; DiscreteEvent (double t, const char event[], double recur_int) : t(t), event(event), recur_int(recur_int) {} DiscreteEvent (double t, const std::string &event, double recur_int) : t(t), event(event), recur_int(recur_int) {} bool operator< (const DiscreteEvent& a) const { return t>a.t; } }; typedef std::priority_queue > deq_type; deq_type deq; public: ///Schedule \a event at \a t void insert(double t, const char event[], double recur_int=0){ assert(recur_int>=0); //assert(empty() || t>=current_time()); deq.push(DiscreteEvent(t,event,recur_int)); } ///Schedule \a event at \a t void insert(double t, const std::string &event, double recur_int=0){ assert(recur_int>=0); //assert(empty() || t>=current_time()); deq.push(DiscreteEvent(t,event,recur_int)); } ///Time the top element of the calendar queue is scheduled for double current_time() const { return deq.top().t; } ///Event indicated by the top element of the calendar queue std::string current_event() const { return deq.top().event; } ///Remove the top element of the calendar queue, and cast it into oblivion void pop() { if(deq.empty()) return; if(deq.top().recur_int>0) insert(current_time()+deq.top().recur_int, deq.top().event, deq.top().recur_int); deq.pop(); } ///Take the top element of the queue and insert a copy of it \a dt ///into the future double reschedule_top(double dt) { deq.push(DiscreteEvent(deq.top().t+dt, deq.top().event, deq.top().recur_int)); return current_time()+dt; } ///Clear all events from the calendar queue void clear() { deq=deq_type(); }; ///Is the calendar queue empty? bool empty() const { return deq.empty(); } }; /** The ArrayState class wraps std::array providing algebra operations used by the Integrator. */ template class ArrayState : public std::array { public: ArrayState() {} ArrayState(const std::vector &init){ for(unsigned int i=0;i::size();++i) std::array::operator[](i)=init[i]; } ArrayState& operator+=(const ArrayState &other) { for(unsigned int i=0;i::size();++i) std::array::operator[](i)+=other[i]; return *this; } }; template ArrayState operator+( const ArrayState &a , double b ){ ArrayState result(a); for(T& i: result) i+=b; return result; } template ArrayState operator+( double a, const ArrayState &b ){ return b+a; } template ArrayState operator+( const ArrayState &a, const ArrayState &b ){ ArrayState result(a); for(unsigned int i=0;i ArrayState operator*( const ArrayState &a, double b ){ ArrayState result(a); for(T& i: result) i*=b; return result; } template ArrayState operator*( double a, const ArrayState &b ){ return b*a; } template ArrayState operator/( const ArrayState &a, double b){ ArrayState result(a); for(T& i: result) i/=b; return result; } template double abs( const ArrayState &a ){ double result=0; for(const T& i: a) result+=abs(i); return result; } template double abs( const ArrayState &a ){ double result=0; for(const double& i: a) result+=std::abs(i); return result; } template std::ostream& operator<<(std::ostream &out, const ArrayState &a){ for(typename ArrayState::const_iterator i=a.begin();i!=a.end();++i) out<<*i<<" "; return out; } /** The Integrator class supplies a simple Euler's method integrator with an adaptive step-size scheme */ template class Integrator { protected: T stateval; typedef std::function dx_type; dx_type dx; double dtval, dtmin, dtmax; double t; int goodsteps; int stepcount; public: Integrator(const T &stateval, dx_type dx, double dtmin, double dtmax); double time() const; T& state(); double dt() const; void dt(double h); void step(); int steps() const; }; template Integrator::Integrator(const T &stateval, dx_type dx, double dtmin, double dtmax) : stateval(stateval), dx(dx), dtmin(dtmin), dtmax(dtmax) { assert(dtmin>0); assert(dtmax>0); assert(dtmin<=dtmax); dtval=dtmin; stepcount=0; t=0; } template double Integrator::time() const { return t; } template T& Integrator::state() { return stateval; } template double Integrator::dt() const { return dtval; } template void Integrator::dt(double h) { assert(h>0); assert(h>=dtmin); assert(h<=dtmax); dtval=h; } template int Integrator::steps() const { return stepcount; } template void Integrator::step() { T e1, e2; ++stepcount; dx(stateval , e1, t); dx(stateval+e1*dtval/2, e2, t+dtval/2.); double abs_e1=abs(stateval+e1*dtval); double abs_e2=abs(stateval+e1*dtval/2+e2*dtval/2); if( (std::abs(abs_e1-abs_e2)/(std::abs(abs_e1+abs_e2)/2))>0.05 ){ stateval+=e1*dtval/2.; t+=dtval/2.; dtval/=2.; if(dtval=8) dtval*=2; if(dtval>dtmax) dtval=dtmax; } /** The EventIntegrator class supplies a simple Euler's method integrator with an adaptive step-size scheme. It also allows for discrete events. It approaches these events with an exponentially-decreasing stepsize and moves away from them in the same fashion. */ template class EventIntegrator : public Integrator { private: CalendarQueue calq; bool at_event; std::string at_event_name; public: EventIntegrator(const T &stateval, typename Integrator::dx_type dx, double dtmin, double dtmax) : Integrator(stateval, dx, dtmin, dtmax) { at_event=false; at_event_name=""; } void insert_event(double t, const std::string &event, double recur_int); void insert_event(double t, const char event[], double recur_int); bool is_event() const; std::string event() const; void step(); }; template void EventIntegrator::insert_event(double t, const std::string &event, double recur_int=0){ calq.insert(t,event,recur_int); } template void EventIntegrator::insert_event(double t, const char event[], double recur_int=0){ calq.insert(t,event,recur_int); } template bool EventIntegrator::is_event() const { return at_event; } template std::string EventIntegrator::event() const { if( is_event() ) return calq.current_event(); else return ""; } template void EventIntegrator::step() { T e1, e2; if(at_event && calq.current_time()==Integrator::t){ at_event_name=calq.current_event(); calq.pop(); return; } at_event=false; ++Integrator::stepcount; if(!calq.empty()){ while( Integrator::dtval > Integrator::dtmin && Integrator::t + Integrator::dtval > calq.current_time() ) { Integrator::dtval/=2.; } if(Integrator::dtval < Integrator::dtmin) Integrator::dtval = Integrator::dtmin; if(Integrator::dtval == Integrator::dtmin && Integrator::t + Integrator::dtval > calq.current_time()){ Integrator::dtval=calq.current_time()-Integrator::t; Integrator::dx(Integrator::stateval, e1, Integrator::t); Integrator::stateval+=e1*Integrator::dtval; Integrator::t=calq.current_time(); Integrator::dtval = Integrator::dtmin; at_event=true; at_event_name=calq.current_event(); calq.pop(); return; } } Integrator::dx(Integrator::stateval , e1, Integrator::t); Integrator::dx(Integrator::stateval + e1*Integrator::dtval/2, e2, Integrator::t+Integrator::dtval/2.); double abs_e1=abs(Integrator::stateval + e1*Integrator::dtval); double abs_e2=abs(Integrator::stateval + e1*Integrator::dtval/2 + e2*Integrator::dtval/2); if( (std::abs(abs_e1-abs_e2)/(std::abs(abs_e1+abs_e2)/2))>0.05 ){ Integrator::stateval+=e1*Integrator::dtval/2.; Integrator::t+=Integrator::dtval/2.; Integrator::dtval/=2.; if(Integrator::dtval::dtmin) Integrator::dtval=Integrator::dtmin; Integrator::goodsteps=0; } else { Integrator::stateval+=e1*Integrator::dtval; Integrator::t+=Integrator::dtval; ++Integrator::goodsteps; } if(Integrator::dtval::dtmax && Integrator::goodsteps>=8) Integrator::dtval*=2; if(Integrator::dtval>Integrator::dtmax) Integrator::dtval=Integrator::dtmax; } } #endif
ċ
simple_integrator-master.zip
(16k)