ni_euler.cc

Go to the documentation of this file.

/////////////////////////////////////////////////////////////////////////////
// ni_euler.cc
//
// Copyright (c) 1991-2004 Petr Peringer
// Copyright (c) 1996-1997 David Leska
//
// This library is licensed under GNU Library GPL. See the file COPYING.
//

//
//  numerical integration: Euler's method
//

////////////////////////////////////////////////////////////////////////////
//  interface
//
#include "simlib.h"
#include "internal.h"
#include "ni_euler.h"
#include <cmath>
#include <cstddef>


////////////////////////////////////////////////////////////////////////////
//  implementation
//

namespace simlib3 {

SIMLIB_IMPLEMENTATION;


////////////////////////////////////////////////////////////////////////////
//  Euler's method
//
/*   Formula:

     y'(t)     = f(t, y);
     y(t+h/2)  = y(t) + h/2 * y'(t);
     y'(t+h/2) = f(t+h/2, y(t+h/2));
     y(t+h)    = y(t+h/2) + h/2 * y'(t+h/2);
     y'(t+h)   = f(t+h, y(t+h));
     err       = h * |y'(t) - y'(t+h/2)|;
*/

void EULER::Integrate(void)
{
  const double err_coef = 0.02; // limits an error range
  static double dthlf;   // half step
  size_t i;   // auxiliary variables for loops to go through list
  Iterator ip, end_it; // of integrators
  static bool DoubleStepFlag; // flag - allow increasing (doubling) the step

  Dprintf((" Euler integration step ")); // print debugging info
  Dprintf((" Time = %g, optimal step = %g", (double)Time, OptStep));

  end_it=LastIntegrator(); // end of container of integrators

  //--------------------------------------------------------------------------
  //  Step of method
  //--------------------------------------------------------------------------

begin_step: // beginning of step

  SIMLIB_StepSize = max(SIMLIB_StepSize, SIMLIB_MinStep); // low step limit

  dthlf = 0.5*SIMLIB_StepSize;     // half step

  SIMLIB_ContractStepFlag = false; // clear reduce step flag
  SIMLIB_ContractStep = 0.5*dthlf; // implicitly reduce to half

  for(ip=FirstIntegrator(),i=0; ip!=end_it; ip++,i++) {
    A[i]   = (*ip)->GetOldDiff();
    (*ip)->SetState((*ip)->GetOldState() + dthlf*(*ip)->GetDiff());   // state y(t+h/2)
  }

  ////////////////////////////////////////////////////////////// 1/2 of step

  _SetTime(Time, SIMLIB_StepStartTime+dthlf);
  SIMLIB_DeltaTime = double(Time)-SIMLIB_StepStartTime;

  SIMLIB_Dynamic();  // compute new state of model                  (1)

  if(StateCond()) { // check on changes of state conditions in 1/2 of step
    goto begin_step;
  }

  bool wasContractStepFlag = SIMLIB_ContractStepFlag; // remember value
  SIMLIB_ContractStepFlag = false; // not reduce step
  SIMLIB_ContractStep = dthlf;     // implicitly reduce to half of step

  StoreState(di, si, xi); // store values in 1/2 of step

  for(ip=FirstIntegrator(),i=0; ip!=end_it; ip++,i++) {
    // difference of differentiations for error estimation
    A[i] -= (*ip)->GetDiff();
    (*ip)->SetState(si[i] + dthlf*(*ip)->GetDiff());
  }

  //////////////////////////////////////////////////////////// end of step

  _SetTime(Time, SIMLIB_StepStartTime + SIMLIB_StepSize);
  SIMLIB_DeltaTime = SIMLIB_StepSize;

  SIMLIB_Dynamic();  // compute new state of model                  (2)

  //--------------------------------------------------------------------------
  //  Check on accuracy of numerical integration, estimate error
  //--------------------------------------------------------------------------

  DoubleStepFlag = true; // allow doubling the step
  SIMLIB_ERRNO = 0; // OK
  for(ip=FirstIntegrator(),i=0; ip!=end_it; ip++,i++) {
    double eerr; // estimated error
    double terr; // greatest allowed error

    eerr = fabs(SIMLIB_StepSize*A[i]); // error estimation
    terr = SIMLIB_AbsoluteError + fabs(SIMLIB_RelativeError*si[i]);

    if(eerr < err_coef*terr) // allowed tolerantion is fulfiled with provision
      continue;

    if(eerr > terr) {        // allowed tolerantion is overfulfiled
      if(SIMLIB_StepSize > SIMLIB_MinStep) {  // reducing step is possible
        SIMLIB_OptStep = 0.5*SIMLIB_StepSize; // halve optimal step
        if(SIMLIB_OptStep < SIMLIB_MinStep) { // limit of optimal step
          SIMLIB_OptStep = SIMLIB_MinStep;
        }
        SIMLIB_StepSize = SIMLIB_OptStep;
        IsEndStepEvent = false;
        goto begin_step; // compute again with smaller step
      }
      // reducing step is unpossible
      SIMLIB_ERRNO++; // requested accuracy cannot be achieved
      _Print("\n Integrator[%lu] ",(unsigned long)i);
      if(SIMLIB_ConditionFlag) // event has been within the step
        break;
    }

    DoubleStepFlag = false;  // disable increasing SIMLIB_OptStep,
                             // accuracy is sufficient, but not well
  } // for

  if(SIMLIB_ERRNO) {
    SIMLIB_warning(AccuracyError);
  }

  //--------------------------------------------------------------------------
  //  Computation is continuing, compute y(t+h)
  //--------------------------------------------------------------------------

  if(wasContractStepFlag) {
    // step reducing has been requested in half step and it is unpossible
    RestoreState(dthlf, di, si, xi); // restore halfstep state
  } else { // go to half step and complete the computation

    GoToState(di, si, xi);

    SIMLIB_StepStartTime += dthlf;
    SIMLIB_DeltaTime = double(Time) - SIMLIB_StepStartTime;

    //-----------------------------------------------------------------------
    //  Analyse system at the end of the step
    //-----------------------------------------------------------------------

    if(StateCond()) { // check on changes of state conditions at end of step
      goto begin_step;
    }
  }

  //--------------------------------------------------------------------------
  //  Results of step have been accepted, take fresh step
  //--------------------------------------------------------------------------

  // increase step, if accuracy was good
  // step increasing is allowed
  // && method is not used to start multi-step method
  if(DoubleStepFlag && !IsStartMode()) {
    SIMLIB_OptStep += SIMLIB_OptStep; // step doubling
  }
  SIMLIB_OptStep = min(SIMLIB_OptStep,SIMLIB_MaxStep); // limit step size

} // EULER::Integrate


}