AnimaL

Tutorial

Documentation

---

Documentation

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

---

explicit_stepper.h

Go to the documentation of this file.

#ifndef ANIMAL_INTEGRATION_EXPLICIT_STEPPER_H
#define ANIMAL_INTEGRATION_EXPLICIT_STEPPER_H

#include <animal/integration/stepper.h>



namespace animal { namespace integration {

// -------------------------------------------------------------
//  
//  Step_Doubling class.
//  
// -------------------------------------------------------------

template <
  class SolverF,
  class HSmallerF,
  class HLargerF,
  class IsDifferentF >

struct Step_Doubling
  : public Stepper_Function<SolverF,HSmallerF,HLargerF,IsDifferentF>
{
  State initial_S_copy;  
  State final_S_eval;    
  
  Derivative initial_D;
  
  unsigned int ngood;  
  unsigned int nbad;   
  
  Step_Doubling()
    {}
  
  Step_Doubling
    (
        const SolverF& slf,
        const HSmallerF& hsf, 
        const HLargerF& hlf,
        const IsDifferentF& isdf
    )
  : Stepper_Function
        <
            SolverF,
            HSmallerF,
            HLargerF,
            IsDifferentF
        >(slf, hsf, hlf, isdf),
    initial_S_copy( slf.curr_size ),
    final_S_eval( slf.curr_size ),
    initial_D( slf.curr_size ),
    ngood(0), nbad(0)
    {}
  
  void resize(const typename Derivative::sizeype size)
    {
      solve.resize(size);
      
      Derivative initial_S_copy(size);
      Derivative final_S_eval(size);
      Derivative initial_D(size);
    }
  
  void operator()(const Model& M,
          State& S,
          const Real t,
          const Real htry, Real& hdid, Real& hnext)
    {
      initial_S_copy = S;
      operator()(M, initial_S_copy, S, t, htry, hdid, hnext);
    }
  void operator()(const Model& M,
          const State& initial_S,
          State& final_S,
          const Real t,
          const Real htry, Real& hdid, Real& hnext)
    {
      bool divide_step = true;
      
      Real h = htry; // Stepsize to be attempted (initial trial value)
      Real h_2;
      
      // Derivative first evaluation
      solve.writeDerivative(M, initial_S, initial_D, t);
      
      while ( divide_step )
    {
      hdid = h; // Stepsize that was actually accomplished
      h_2 = 0.5*h;
      
      // Take a full step
      solve(M, initial_S, final_S_eval, initial_D, t, h);
      
        // Take two half steps
        solve(M, initial_S, final_S, initial_D, t, h_2);
        solve(M, final_S, (t + h_2), h_2);

        if ( isDifferent(M, final_S_eval, final_S, t) )
      {
          h = hsmaller(M, initial_S, t, h);
#if DEBUG
          nbad++;
#endif
      }
      else
        {
          divide_step = false;
          h = hlarger(M, initial_S, t, h);
#if DEBUG
          ngood++;
#endif
        }
    }
      
      hnext = h; // Estimated next stepsize
    }
  
}; // struct Step_Doubling



// -------------------------------------------------------------
//  
//  Step_Back_And_Forth class.
//  
// -------------------------------------------------------------

template <
  class SolverF,
  class HSmallerF,
  class HLargerF,
  class IsDifferentF >

struct Step_Back_And_Forth
  : public Step_Doubling<SolverF,HSmallerF,HLargerF,IsDifferentF>
{
  Step_Back_And_Forth(const SolverF& slf,
              const HSmallerF& hsf, const HLargerF& hlf,
              const IsDifferentF& isdf)
  : Step_Doubling<SolverF,HSmallerF,HLargerF,IsDifferentF>(slf, hsf, hlf, isdf)
    {}
  
  void operator()(const Model& M,
          State& S,
          const Real t,
          const Real htry, Real& hdid, Real& hnext)
    {
      initial_S_copy = S;
      operator()(M, initial_S_copy, S, t, htry, hdid, hnext);
    }
  void operator()(const Model& M,
          const State& initial_S,
          State& final_S,
          const Real t,
          const Real htry, Real& hdid, Real& hnext)
    {
      bool divide_step = true;
      
      Real h = htry; // Stepsize to be attempted (initial trial value)
      
      // Derivative first evaluation
      solve.writeDerivative(M, initial_S, initial_D, t);
      
      while ( divide_step )
    {
      hdid = h; // Stepsize that was actually accomplished
          
          // Take a step
          solve(M, initial_S, final_S, initial_D, t, h);
          
          if ( isDifferent(M, initial_S, final_S, t) )
        {
          h = hsmaller(M, initial_S, t, h);
#if DEBUG
          nbad++;
#endif
        }
      else
        {
          divide_step = false;
          h = hlarger(M, initial_S, t, h);
#if DEBUG
          ngood++;
#endif
        }
    }
      
      hnext = h; // Estimated next stepsize
    }
  
}; // struct Step_Back_And_Forth

} } // namespace animal { namespace integration {



#endif // ANIMAL_INTEGRATION_EXPLICIT_STEPPER_H

---