AnimaL

Tutorial

Documentation

---

Documentation

• Overview
• Modules
• Namespace List
• Class Hierarchy
• Alphabetical List
• Compound List
• File List
• Namespace Members
• Compound Members
• File Members
• Examples

---

integrator.h

Go to the documentation of this file.

#ifndef ANIMAL_INTEGRATOR_HPP
#define ANIMAL_INTEGRATOR_HPP

#include <animal/linear.h>
#include <animal/vectorVec.h>


namespace animal {

template<class FS, class IS, class D, class Real>
inline void v_eq_euler_step( FS& final, const IS& initial, Real h, const D& derivative )
{
    //v_eq(final,initial,h,derivative);
    v_eq(final.positions, initial.positions, h, derivative.velocities );
    v_eq(final.velocities, initial.velocities, h, derivative.accelerations );
}

template<class State, class Deriv, class Derivator, class Real>
inline void integrate_euler( State& s, Real h, Derivator* eq_deriv, Real t, Deriv& d )
{
    (*eq_deriv)(d,s,t);
    v_eq_euler_step(s,s,h,d);
}

template<class State, class Deriv, class Derivator, class StateAux, class Real>
inline void integrate_rk2( State& s, Real h, Derivator* eq_deriv, Real t, Deriv& d, StateAux& s1 )
{
    (*eq_deriv)(d,s,t);
    v_eq_euler_step(s1,s,h/2,d);
    (*eq_deriv)(d,s1,t+(h/2));
    v_eq_euler_step(s,s,h,d);
}


template<class State, class Deriv, class Derivator, class StateAux, class Real>
inline void integrate_rk4( State& s, Real h, Derivator* eq_deriv, Real t, Deriv& d1, Deriv& d2, Deriv& d3, Deriv& d4, StateAux& s1 )
{
    (*eq_deriv)(d1,s,t);
    v_eq_euler_step(s1,s,h/2,d1);
    (*eq_deriv)(d2,s1,t+(h/2));
    v_eq_euler_step(s1,s,h/2,d2);
    (*eq_deriv)(d3,s1,t+(h/2));
    v_eq_euler_step(s1,s,h,d3);
    (*eq_deriv)(d4,s1,t+h);
    
    // accumulate composed derivative in d2
    v_peq(d2,d3); // d2+d3
    v_teq(d2,2.);  // 2(d2+d3)
    v_peq(d2,d1); // d1+2(d2+d3)
    v_peq(d2,d4); // d1+2(d2+d3)+d4
    
    v_eq_euler_step(s,s,h/6,d2);
}


template<class State, class Deriv, class Derivator, class StateAux, class Real>
inline void integrate_modmid( State& s, Real h, Derivator *eq_deriv, Real t, int n, Deriv& d1, StateAux& s1, StateAux& s2, StateAux& s3 )
{
    Real hsub = h/n;       // h: total step to be made
    Real hhsub = 2.0*hsub; // hsub: substep value

    v_eq(s1,s);

    (*eq_deriv)(d1,s,t);
    v_eq_euler_step(s2,s,hsub,d1); // First substep taken
    
    // take remaining substeps minus one
    int i;
    for( i=0; i<n-1; ++i )
    {
        (*eq_deriv)(d1,s2,t+i*hsub);
        v_eq_euler_step(s3,s1,hhsub,d1);
        v_eq(s1,s2);
        v_eq(s2,s3);
    }
    
    //  Take last substep
    (*eq_deriv)(d1,s2,t+i*hsub);
    v_eq_euler_step(s3,s1,hhsub,d1);
    
    // Take another step
    (*eq_deriv)(d1,s3,t+h);
    v_eq_euler_step(s3,s3,hsub,d1);
    
    // result = average(s2,s3)
    v_peq(s2,s3);
    Real unDemi = 0.5;
    v_eq_ab( s,unDemi,s2 );
}



template<class State, class Deriv, class Derivator, class Real>
class ExplicitIntegrator
{
    Derivator* eq_deriv;  
public:
    ExplicitIntegrator( Derivator* d=0 )
        : eq_deriv(d)
    {}
    
    template<class States>
    void euler(States& s, Real h, Real t=0)
    {
        assert(eq_deriv);
        d1.resize( size(s) );
        integrate_euler(s,h,eq_deriv,t,d1);
    }

    template<class States>
    void rk2(States& s, Real h, Real t=0)
    {
        assert(eq_deriv);
        d1.resize( size(s) );
        s1.resize( size(s) );
        integrate_rk2(s,h,eq_deriv,t,d1,s1);
    }

    template<class States>
    void rk4(States& s, Real h, Real t=0)
    {
        assert(eq_deriv);
        d1.resize( size(s) );
        d2.resize( size(s) );
        d3.resize( size(s) );
        d4.resize( size(s) );
        s1.resize( size(s) );
        integrate_rk4(s,h,eq_deriv,t,d1,d2,d3,d4,s1);
    }

    template<class States>
    void modmid(States& s, Real h, int n, Real t=0)
    {
        assert(eq_deriv);
        d1.resize( size(s) );
        s1.resize( size(s) );
        s2.resize( size(s) );
        s3.resize( size(s) );
        integrate_modmid(s,h,eq_deriv,t,n,d1,s1,s2,s3);
    }

protected:
    animal::vector<Deriv> d1;  
    animal::vector<Deriv> d2;  
    animal::vector<Deriv> d3;  
    animal::vector<Deriv> d4;  
    animal::vector<State> s1;  
    animal::vector<State> s2;  
    animal::vector<State> s3;  
};





//==========================================================================
// Special case of ExplicitIntegrator for second-degree systems
// e.g. (position, velocity) -> (velocity,acceleration)

template<class Poss, class Vels>
struct States {
    Poss& positions;  
    Vels& velocities;  

    States( Poss& pos, Vels& vel )
        : positions(pos)
        , velocities(vel)
    {}

    void resize( int s ){
        velocities.resize(s);
        positions.resize(s);
    }
};

template<class P,class V> inline 
void v_eq(States<P,V>& s, const States<P,V>& t ){
    v_eq(s.positions,t.positions);
    v_eq(s.velocities,t.velocities);
}

template<class P,class V> inline 
void v_peq(States<P,V>& s, const States<P,V>& t ){
    v_peq(s.positions,t.positions);
    v_peq(s.velocities,t.velocities);
}

template<class P,class V,class Real> inline 
void v_eq_ab(States<P,V>& s, const Real& a, const States<P,V>& b ){
    v_eq_ab(s.positions,a,b.positions);
    v_eq_ab(s.velocities,a,b.velocities);
}

template<class Vels,class Accs>
struct Derivs {
    Accs& accelerations;   
    Vels& velocities;           

    Derivs( Vels& vel, Accs& acc )
        : accelerations(acc)
        , velocities(vel)
    {}

    void resize( int s ){
        velocities.resize(s);
        accelerations.resize(s);
    }
};

template<class V,class A> inline 
void v_peq( Derivs<V,A>& d, const Derivs<V,A>& e ){
    v_peq(d.velocities,e.velocities);
    v_peq(d.accelerations,e.accelerations);
}

template<class V,class A,class Real> inline 
void v_teq( Derivs<V,A>& d, const Real& e ){
    v_teq(d.velocities,e);
    v_teq(d.accelerations,e);
}


template<class Position, class Velocity, class ComputeAcceleration, class Real, class Acceleration=Velocity>
class ExplicitIntegrator2
{
    typedef animal::VectorVec<Position> Positions;
    
    typedef animal::VectorVec<Velocity> Velocities;
    
    typedef animal::VectorVec<Acceleration> Accelerations;
    
    
    struct Derivator {
        ComputeAcceleration& computeAcceleration;
        
        Derivator(ComputeAcceleration& compAcc)
            : computeAcceleration(compAcc)
        {}
        
        template<class P,class V,class A>
        void operator() (Derivs<V,A>& deriv, const States<P,V>& state, Real /*time not used*/)
        {
            computeAcceleration( deriv.accelerations, state.positions, state.velocities );
            v_eq( deriv.velocities, state.velocities );
        }
    };
    
    Derivator eq_deriv;  
public:
    ExplicitIntegrator2( ComputeAcceleration& d )
        : eq_deriv(d)
        , d1(dv1,da1)
        , d2(dv2,da2)
        , d3(dv3,da3)
        , d4(dv4,da4)
        , s1(sp1,sv1)
        , s2(sp2,sv2)
        , s3(sp3,sv3)
    {}
    
    template< class P, class V >
    void euler(P& pos, V& vel, Real h, Real t=0)
    {
        d1.resize( size(pos) );
        States<P,V> states(pos,vel);
        //cerr<<"integrator2::euler: states.pos= "<< v_output(states.positions) <<endl;
        //cerr<<"integrator2::euler: states.vel= "<< v_output(states.velocities) <<endl;
        integrate_euler(states,h,&eq_deriv,t,d1);
    }

    template<class P, class V>
    void rk2(P& pos, V& vel, Real h, Real t=0)
    {
        d1.resize( size(pos) );
        s1.resize( size(pos) );
        States<P,V> states(pos,vel);
        integrate_rk2(states,h,&eq_deriv,t,d1,s1);
    }

    template<class P, class V>
    void rk4(P& pos, V& vel, Real h, Real t=0)
    {
        d1.resize( size(pos) );
        d2.resize( size(pos) );
        d3.resize( size(pos) );
        d4.resize( size(pos) );
        s1.resize( size(pos) );
        States<P,V> states(pos,vel);
        integrate_rk4(states,h,&eq_deriv,t,d1,d2,d3,d4,s1);
    }

    template<class P, class V>
    void modmid(P& pos, V& vel, Real h, int n, Real t=0)
    {
        d1.resize( size(pos) );
        s1.resize( size(pos) );
        s2.resize( size(pos) );
        s3.resize( size(pos) );
        States<P,V> states(pos,vel);
        integrate_modmid(states,h,&eq_deriv,t,n,d1,s1,s2,s3);
    }

protected:
    
    typedef Derivs<Velocities,Accelerations> Der;  
    
    Velocities dv1;  
    Accelerations da1;  
    Der d1; 
    
    Velocities dv2;  
    Accelerations da2;  
    Der d2; 
    
    Velocities dv3;  
    Accelerations da3;  
    Der d3; 
    
    Velocities dv4;  
    Accelerations da4;  
    Der d4; 
    
    
    typedef States<Positions,Velocities> Sta;  
    
    Positions sp1;   
    Velocities sv1;   
    Sta s1;        

    Positions sp2;   
    Velocities sv2;   
    Sta s2;        

    Positions sp3;   
    Velocities sv3;   
    Sta s3;        

};


} // end namespace animal



#endif 

---