ODESystem.h

Go to the documentation of this file.

#ifndef ODESYSTEM_H_
#define ODESYSTEM_H_

#include <gsl/gsl_errno.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_odeiv.h>

#include <string>
using std::string;

#include <iostream>
using std::cerr;
using std::endl;

namespace odeiv {

    class ODESystem {

        public:
            ODESystem( size_t dim ) {
                gsl_sys.function = __int__dydt__;
                gsl_sys.jacobian = NULL;
                gsl_sys.dimension = dim;
                gsl_sys.params = (void *)this;
            };

            void setSystem( ODESystem *sys ) {
                gsl_sys.params = (void *)sys;
            };

            virtual ~ODESystem() { }
            ;

            virtual void derivatives( double t, const double y[], double f[] ) = 0;

            int dimension() const {
                return gsl_sys.dimension;
            };

        private:
            friend class GslStepper;
            friend class AdaptiveStepSolver;
            gsl_odeiv_system gsl_sys;
            static int __int__dydt__( double t, const double y[], double f[], void *params ) {
                (( ODESystem* )( params ))->derivatives( t, y, f );
                return GSL_SUCCESS;
            }
    };


    class Stepper {

        public:

            Stepper( ) {
                /* NOOP */
            };

            virtual ~Stepper( ) {
                /* NOOP */
            };

            //        virtual const gsl_odeiv_step_type *type()=0;
            virtual void reset( ) = 0;
            virtual string name() = 0;
            virtual int order() = 0;
            virtual void apply( double t, double h, double y[], double yerr[], const double dydt_in[], double dydt_out[], ODESystem  & sys ) = 0;

    };

    class GslStepperType {
        public:
            virtual ~GslStepperType() {}
            ;
            virtual const gsl_odeiv_step_type *gsl_type() const =0;
    };

    class RungeKutta2 : public GslStepperType {
        public:
            virtual const gsl_odeiv_step_type *gsl_type() const {
                return gsl_odeiv_step_rk2;
            };
    };

    class RungeKutta4 : public GslStepperType {
        public:
            virtual const gsl_odeiv_step_type *gsl_type() const  {
                return gsl_odeiv_step_rk4;
            };
    };

    class GslStepper : public Stepper {

        public:
            GslStepper( ODESystem & sys, GslStepperType const& step ) : int_step( NULL ) {
                int_step = gsl_odeiv_step_alloc( step.gsl_type(), sys.dimension() );
            };

            virtual ~GslStepper( ) {
                if( int_step )
                    gsl_odeiv_step_free( int_step );
            };

            virtual void reset( ) {
                gsl_odeiv_step_reset( int_step );
            };

            virtual string name() {
                return string( gsl_odeiv_step_name( int_step ) );
            };

            virtual int order() {
                return gsl_odeiv_step_order( int_step );
            };

            virtual void apply( double t, double h, double y[], double yerr[], const double dydt_in[], double dydt_out[], ODESystem  & sys ) {
                gsl_odeiv_step_apply( int_step, t, h, y, yerr, dydt_in, dydt_out, &sys.gsl_sys );
            };

        protected:
            friend class AdaptiveStepSolver;
            gsl_odeiv_step *int_step;

    };


    class Solver {

        public:

            Solver( ODESystem & sys, GslStepperType const& step ) : stepper( sys, step ) {
                _y_ = new double[sys.dimension()];
            };

            virtual ~Solver() {
                if ( _y_ )
                    delete _y_;
            };

            virtual void reset( double dt, ODESystem  & sys, const double y0[] ) = 0;

            virtual void advance( double t, double dt, ODESystem  & sys ) = 0;

            double y( size_t i ) {
                return _y_[ i ];
            };

        protected:
            GslStepper stepper;
            double *_y_;
    };

    class FixedStepSolver : public Solver {

        public:
            FixedStepSolver( ODESystem & sys, GslStepperType const& step );
            virtual ~FixedStepSolver();
            //virtual void init( ODESystem &sys );
            virtual void reset( double dt, ODESystem & sys, const double y0[] );
            virtual void advance( double t, double dt, ODESystem  & sys );

        private:
            double *dydt_in, *dydt_out, *y_err;
    };


    class AdaptiveStepSolver : public Solver {

        public:
            AdaptiveStepSolver( ODESystem & sys, GslStepperType const& step );
            virtual ~AdaptiveStepSolver();
            // virtual void init( ODESystem & sys );
            virtual void reset( double dt, ODESystem & sys, const double y0[] );
            virtual void advance( double t, double dt, ODESystem  & sys );

        private:
            gsl_odeiv_control *int_control;
            gsl_odeiv_evolve *int_evolve;
            double h;
    };


}
#endif /*ODESYSTEM_H_*/

---