ODESystem.cpp

Go to the documentation of this file.

#include "ODESystem.h"
using namespace odeiv;

#include "PCSIMException.h"

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

FixedStepSolver::FixedStepSolver( ODESystem & sys, GslStepperType const& step ) : Solver( sys, step ) {
    dydt_in = new double[ sys.dimension() ];
    dydt_out = new double[ sys.dimension() ];
    y_err = new double[ sys.dimension() ];
}

FixedStepSolver::~FixedStepSolver() {
    if ( dydt_in )
        delete dydt_in;
    if ( dydt_out )
        delete dydt_out;
    if ( y_err )
        delete y_err;
}

void FixedStepSolver::reset( double dt, ODESystem & sys, const double y0[] ) {
    for ( int i = 0; i < sys.dimension(); i++ ) {
        _y_[ i ] = y0[ i ];
    }
    sys.derivatives( 0.0, _y_, dydt_in );
    stepper.reset( );
}

void FixedStepSolver::advance( double t, double dt, ODESystem & sys ) {
    // cerr << "FIXED STEP: " << dt ;
    stepper.apply( t, dt, _y_, y_err, dydt_in, dydt_out, sys );
    for ( int i = 0; i < sys.dimension(); i++ )
        dydt_in[ i ] = dydt_out[ i ];
    // cerr << " : Vm = " << _y_[ 0 ] << endl;
}


AdaptiveStepSolver::AdaptiveStepSolver( ODESystem & sys, GslStepperType const& step ) : Solver( sys, step ) {
    int_control = gsl_odeiv_control_y_new( 1e-6, 1e-6 );
    int_evolve = gsl_odeiv_evolve_alloc( sys.dimension() );
    h = 1e-6;
    //cerr << "AdaptiveStepSolver::AdaptiveStepSolver" << endl;
}


AdaptiveStepSolver::~AdaptiveStepSolver() {
  //cerr << "DESTRUCTOR AdaptiveStepSolver::AdaptiveStepSolver" << endl;
    if ( int_control != NULL )
        gsl_odeiv_control_free( int_control );
    if ( int_evolve != NULL )
        gsl_odeiv_evolve_free( int_evolve );
}

void AdaptiveStepSolver::reset( double dt, ODESystem & sys, const double y0[] ) {
    // cerr << "ADAPT RESET: " << dt << " ep=" << int_evolve << " sp=" << stepper->int_step << "sd=" << ( stepper->int_step ) ->dimension << endl;
    for ( int i = 0; i < sys.dimension(); i++ ) {
        _y_[ i ] = y0[ i ];
    }
    stepper.reset( );
    gsl_odeiv_evolve_reset( int_evolve );
    h = dt / 100.0;
    // cerr << "ADAPT RESET: " << dt << " ep=" << int_evolve << endl;
}

void AdaptiveStepSolver::advance( double t, double dt, ODESystem & sys ) {
    // cerr << "ADAPT STEP: " << " t=" << t << " dt=" << dt << endl ;
    double t1 = t + dt;
    //cerr << "ADAPT STEP: dt=" << dt << ", h=" << h << ", y = ";
    //for ( int i = 0; i < sys.dimension(); i++ ) {
    //  cerr << _y_[ i ] << ", " ;
    //}
    while ( t < t1 ) { 
        // cerr << "while ( t < t1 ) { " << dt << " ep=" << int_evolve->dimension << " sp=" << ( stepper->int_step ) ->dimension << endl;
        if ( gsl_odeiv_evolve_apply( int_evolve, int_control, stepper.int_step, &sys.gsl_sys, &t, t1, &h, _y_ ) != GSL_SUCCESS ) {
            throw( PCSIM::Exception( "AdaptiveStepSolver::advance", "ODE Solver error!" ) );
        }
        // cerr << " -> " << t << " ";
    }
    //cerr << ": Vm(" << t1 << ") = " << _y_[ 0 ] << endl;
}

---