AnimaL

Tutorial

Documentation

---

Documentation

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

---

explicit_solver.h

Go to the documentation of this file.

#ifndef ANIMAL_INTEGRATION_EXPLICIT_SOLVER_H
#define ANIMAL_INTEGRATION_EXPLICIT_SOLVER_H

#include <animal/integration/solver.h>



namespace animal { namespace integration {

// ------------------------------------------------------------
//  
//  Euler class.
//  
// ------------------------------------------------------------

template 
<
    class DerivativeF,
    class StepF,
    class TraitsT
>

struct Euler : public Solver< DerivativeF, StepF, TraitsT >
{
  
//  typedef DerivativeF Derivative;   ///< derivative functor
//  typedef StepF Step;               ///< basic step functor
//  typedef TraitsT Traits;           ///< traits

  typedef TraitsT                     Traits;

  typedef typename Traits::Real       Real;

  typedef typename Traits::State      State;

  typedef typename Traits::CopyState      CopyState;

  typedef typename Traits::Derivative Derivative;

  typedef typename Traits::SizeType   SizeType;


    

  Euler
    (
        const DerivativeF& df = DerivativeF(), 
        const StepF& sf = StepF(),
        const TraitsT& tr = Traits()
    )
  : Solver< DerivativeF, StepF, TraitsT >( df, sf, tr )
    {}
    
    virtual ~Euler(){}
    
    


    void doStep
    (
          const State& initial_S,
          State& final_S,
          const Derivative& initial_D,
          const Real, // t
            const Real h
    )
  {
    applyStep( initial_S, final_S, initial_D, h );
  }
    
  

protected:


  void setSize( const SizeType ){}
    

}; // struct Euler



// ------------------------------------------------------
//  
//  Runge_Kutta_2 class.
//  
// ------------------------------------------------------

template 
<
  class DerivativeF,
  class StepF,
  class TraitsT
>

struct Runge_Kutta_2 : public Solver< DerivativeF,StepF,TraitsT >
{

  typedef TraitsT                     Traits;

  typedef typename Traits::Real       Real;

  typedef typename Traits::State      State;

  typedef typename Traits::CopyState      CopyState;

  typedef typename Traits::Derivative Derivative;

  typedef typename Traits::SizeType   SizeType;



  Runge_Kutta_2
    (
        const DerivativeF& df = DerivativeF(), 
        const StepF& sf = StepF(),
        const Traits& tr = Traits()
    )
  : Solver< DerivativeF,StepF,TraitsT >( df, sf, tr )
  {}
    
    virtual ~Runge_Kutta_2(){}
    
  


  void doStep
    (
        const State& initial_S,
        State& final_S,
        const Derivative& D,
        const Real t, 
        const Real h
    )
  {
    Real h_2 = 0.5*h;

    applyStep(initial_S, tmp_S, D, h_2);

    writeDerivative( tmp_S, D1, (t + h_2) );
    applyStep( initial_S, final_S, D1, h );
  }
    
  
protected:


  CopyState tmp_S;
  
  void setSize( const SizeType size )
  {
        traits.resize( tmp_S, size );
  }
    
  
}; // class Runge_Kutta_2


//  /** Apply a second-order Runge-Kutta integration step
//      \param finalState The final state
//      \param initialState The initial state
//      \param doDerivate A functor which computes a derivative given a state
//      \param h The time step
//      \param derivative An auxiliary value to store a derivative
//      \param doStep A functor which applies a basic Euler step
//  */
//  template< class FS, class IS, class DR, class DV, class S>
//  void rk2Step( FS& finalState, const IS& initialState, DR& doDerivate, double h, DV& derivative, S& doStep )
//  {
//      doDerivate( derivative, initialState );
//      doStep( finalState, initialState, derivative, h/2 );  // half step
//      doDerivate( derivative, finalState );                 // derivative there
//      doStep( finalState, initialState, derivative, h );
//  }



// -------------------------------------------------------
//  
//  Runge_Kutta_4 class.
//  
// -------------------------------------------------------

template 
<
  class DerivativeF,
  class StepF,
  class TraitsT
>

struct Runge_Kutta_4 : public Solver< DerivativeF,StepF,TraitsT >
{
  
  typedef TraitsT                     Traits;

  typedef typename Traits::Real       Real;

  typedef typename Traits::State      State;

  typedef typename Traits::CopyState      CopyState;

  typedef typename Traits::Derivative Derivative;

  typedef typename Traits::SizeType   SizeType;



  Runge_Kutta_4
    (
        const DerivativeF& df = DerivativeF(), 
        const StepF& sf = StepF(),
        const Traits& tr = Traits()
    )
  : Solver< DerivativeF, StepF, TraitsT >( df, sf, tr )
  {}
    
    virtual ~Runge_Kutta_4(){}
    
  

  

  void doStep
    (
        const State& initial_S,
        State& final_S,
        const Derivative& initial_D,
        const Real t, 
        const Real h
    )
  {
        Real h_2 = 0.5*h;

        applyStep(initial_S, tmp_S, initial_D, h_2);

        writeDerivative ( tmp_S, D2, (t + h_2) );
        applyStep( initial_S, tmp_S, D2, h_2 );

        writeDerivative ( tmp_S, D3, (t + h_2) );
        applyStep( initial_S, tmp_S, D3, h );

        writeDerivative ( tmp_S, D4, (t + h) );
        
        // weighted sum of derivatives: final_D = initial_D + 2*(D2+d3) + D4
        traits.rk4Step( initial_D, D2, D3, D4, final_D );

        applyStep(initial_S, final_S, final_D, h/6 );

  }

  
protected:


  CopyState tmp_S;
  
  Derivative D2, D3, D4;
  
  Derivative final_D;

  void setSize(const SizeType size)
  {
        traits.resize( tmp_S,  size );
        traits.resize( D2, size );
        traits.resize( D3, size );
        traits.resize( D4, size );
        traits.resize( final_D, size );
  }
    

}; // class Runge_Kutta_4



// --------------------------------------
//  
//  Modified_Midpoint class.
//  
// --------------------------------------

template 
<
  int N,
  class DerivativeF,
  class StepF,
  class TraitsT
>

struct Modified_Midpoint : public Solver< DerivativeF, StepF, TraitsT >
{
  
  
  typedef TraitsT                     Traits;

  typedef typename Traits::Real       Real;

  typedef typename Traits::State      State;

  typedef typename Traits::CopyState      CopyState;

  typedef typename Traits::Derivative Derivative;

  typedef typename Traits::SizeType   SizeType;



  Modified_Midpoint
    (
        const DerivativeF& df = DerivativeF(), 
        const StepF& sf = StepF(),
        const Traits& tr = Traits()
    )
  : Solver< DerivativeF, StepF, TraitsT >( df, sf, tr )
  {}
    
  

  void doStep
    (
        const State& initial_S,
        State& final_S,
        const Derivative& initial_D,
        const Real t, 
        const Real h
    )
  {

    Real hsub = h/N;       // h: total step to be made
    Real hhsub = 2.0*hsub; // hsub: substep value

    traits.copy( initial_S, S1 );

    applyStep( S1, S2, initial_D, hsub ); // First substep taken

    // Take remaining substeps minus one
    for( int m = 1; m < (N - 1); m++ )
        {
            writeDerivative( S2, D1, (t + m*hsub) );
            applyStep(S1, S3, D1, hhsub);
            S1 = S2;
            S2 = S3;
        }

    // Take last substep
    writeDerivative( S2, D1, (t + (N - 1)*hsub) );
    applyStep( S1, S3, D1, hhsub );

    // Take full step
    writeDerivative( S3, D1, (t + h) );
    applyStep( S3, S3, D1, hsub );
        traits.average( S2, S3, final_S );
        
//          cout<<" performed modified_midpoint with N = " << N << endl;

  }
    

  
protected:


  CopyState S1, S2, S3;
  
  void setSize( const SizeType size )
  {
        traits.resize( S1, size );
        traits.resize( S2, size );
        traits.resize( S3, size );
  }
    
  
}; // class Modified_Midpoint



// --------------------------------------
//  
//  Modified_Midpoint_V class.
//  
// --------------------------------------

template 
<
  class NStepF,
  class DerivativeF,
  class StepF,
  class TraitsT
>

struct Modified_Midpoint_V : public Solver< DerivativeF, StepF, TraitsT >
{
  
  typedef TraitsT                     Traits;

  typedef typename Traits::Real       Real;

  typedef typename Traits::State      State;

  typedef typename Traits::CopyState      CopyState;

  typedef typename Traits::Derivative Derivative;

  typedef typename Traits::SizeType   SizeType;


  NStepF nstep;  
  

  Modified_Midpoint_V
    (
        const NStepF& nsf = NStepF(),
        const DerivativeF& df = DerivativeF(), 
        const StepF& sf = StepF(),
        const Traits& tr = Traits()
    )
  : Solver< DerivativeF, StepF, TraitsT >( df, sf, tr )
    , nstep(nsf)
  {}
    
    virtual ~Modified_Midpoint_V(){}
    
  

  void doStep
    (
        const State& initial_S,
        State& final_S,
        const Derivative& initial_D,
        const Real t, 
        const Real h
    )
  {
    unsigned int N = nstep( initial_S, t );

    Real hsub = h/N;       // h: total step to be made
    Real hhsub = 2.0*hsub; // hsub: substep value

    S1 = initial_S;

    applyStep( S1, S2, initial_D, hsub ); // First substep taken

    // Take remaining substeps minus one
    for( unsigned int m = 1; m < (N - 1); m++ )
        {
            writeDerivative( S2, D1, (t + m*hsub) );
            applyStep(S1, S3, D1, hhsub);
            S1 = S2;
            S2 = S3;
        }

    // Take last substep
    writeDerivative( S2, D1, (t + (N - 1)*hsub) );
    applyStep( S1, S3, D1, hhsub );

    // Take full step
    writeDerivative( S3, D1, (t + h) );
    applyStep( S3, S3, D1, hsub );
        traits.average( S2, S3, final_S );
        
//          cout<<" performed modified_midpoint with N = " << N << endl;

  }
    

  
protected:


  State S1, S2, S3;
  
  void setSize( const SizeType size )
  {
        traits.resize( S1, size );
        traits.resize( S2, size );
        traits.resize( S3, size );
  }
    
  
}; // class Modified_Midpoint






template <class TraitsT>

struct NStep_Function
{
  typedef typename TraitsT::Real     Real;     
  typedef typename TraitsT::State    State;    
  typedef typename TraitsT::Numerics Numerics; 
  
  unsigned int operator()
    (
        const State& S,
        const Real t
    ) 
    const;
    
};

} // end namespace integration
} // end namespace animal



#endif // ANIMAL_INTEGRATION_EXPLICIT_SOLVER_H

---