rkf45.c

Runge-Kutta-Fehlberg method

//  Advance an approximate solution over one step of length H.
//
      fehl_d ( f, neqn, y, *t, h, yp, f1, f2, f3, f4, f5, f1 );
      nfe = nfe + 5;
//
//  Compute and test allowable tolerances versus local error estimates
//  and remove scaling of tolerances.  The relative error is
//  measured with respect to the average of the magnitudes of the
//  solution at the beginning and end of the step.
//
      eeoet = 0.0E+00;
 
      for ( k = 0; k < neqn; k++ )
      {
        et = fabs ( y[k] ) + fabs ( f1[k] ) + ae;

        if ( et <= 0.0E+00 )
        {
          delete [] f1;
          delete [] f2;
          delete [] f3;
          delete [] f4;
          delete [] f5;
          return 5;
        }

        ee = fabs 
        ( ( -2090.0E+00 * yp[k] 
          + ( 21970.0E+00 * f3[k] - 15048.0E+00 * f4[k] ) 
          ) 
        + ( 22528.0E+00 * f2[k] - 27360.0E+00 * f5[k] ) 
        );

        eeoet = d_max ( eeoet, ee / et );

      }

      esttol = fabs ( h ) * eeoet * scale / 752400.0E+00;

      if ( esttol <= 1.0E+00 )
      {
        break;
      }
//
//  Unsuccessful step.  Reduce the stepsize, try again.
//  The decrease is limited to a factor of 1/10.
//
      hfaild = true;
      output = false;

      if ( esttol < 59049.0E+00 )
      {
        s = 0.9E+00 / pow ( esttol, 0.2E+00 );
      }
      else
      {
        s = 0.1E+00;
      }

      h = s * h;

      if ( fabs ( h ) < hmin )
      {
        kflag = 6;
        delete [] f1;
        delete [] f2;
        delete [] f3;
        delete [] f4;
        delete [] f5;
        return 6;
      }

    }
//
//  We exited the loop because we took a successful step.  
//  Store the solution for T+H, and evaluate the derivative there.
//
    *t = *t + h;
    for ( i = 0; i < neqn; i++ )
    {
      y[i] = f1[i];
    }
    f ( *t, y, yp );
    nfe = nfe + 1;
//
//  Choose the next stepsize.  The increase is limited to a factor of 5.
//  If the step failed, the next stepsize is not allowed to increase.
//
    if ( 0.0001889568E+00 < esttol )
    {
      s = 0.9E+00 / pow ( esttol, 0.2E+00 );
    }
    else
    {
      s = 5.0E+00;
    }

    if ( hfaild )
    {
      s = d_min ( s, 1.0E+00 );
    }

    h = d_sign ( h ) * d_max ( s * fabs ( h ), hmin );
//
//  End of core integrator
//
//  Should we take another step?
//
    if ( output )
    {
      *t = tout;
      delete [] f1;
      delete [] f2;
      delete [] f3;
      delete [] f4;
      delete [] f5;
      return 2;
    }

    if ( flag <= 0 )
    {
      delete [] f1;
      delete [] f2;
      delete [] f3;
      delete [] f4;
      delete [] f5;
      return (-2);
    }

  }
# undef MAXNFE
}
//******************************************************************************

int rkf45_s ( void f ( float t, float y[], float yp[] ), int neqn, 
  float y[], float yp[], float *t, float tout, float *relerr, float abserr, 
  int flag )

//******************************************************************************
//
//  Purpose:
//
//    RKF45_S carries out the Runge-Kutta-Fehlberg method (single precision).
//
//  Discussion:
//
//    This routine is primarily designed to solve non-stiff and mildly stiff
//    differential equations when derivative evaluations are inexpensive.
//    It should generally not be used when the user is demanding
//    high accuracy.
//
//    This routine integrates a system of NEQN first-order ordinary differential
//    equations of the form:
//
//      dY(i)/dT = F(T,Y(1),Y(2),...,Y(NEQN))
//
//    where the Y(1:NEQN) are given at T.
//
//    Typically the subroutine is used to integrate from T to TOUT but it
//    can be used as a one-step integrator to advance the solution a
//    single step in the direction of TOUT.  On return, the parameters in
//    the call list are set for continuing the integration.  The user has
//    only to call again (and perhaps define a new value for TOUT).
//
//    Before the first call, the user must 
//
//    * supply the subroutine F(T,Y,YP) to evaluate the right hand side;
//      and declare F in an EXTERNAL statement;
//
//    * initialize the parameters:
//      NEQN, Y(1:NEQN), T, TOUT, RELERR, ABSERR, FLAG.
//      In particular, T should initially be the starting point for integration,
//      Y should be the value of the initial conditions, and FLAG should 
//      normally be +1.
//
//    Normally, the user only sets the value of FLAG before the first call, and
//    thereafter, the program manages the value.  On the first call, FLAG should
//    normally be +1 (or -1 for single step mode.)  On normal return, FLAG will
//    have been reset by the program to the value of 2 (or -2 in single 
//    step mode), and the user can continue to call the routine with that 
//    value of FLAG.
//
//    (When the input magnitude of FLAG is 1, this indicates to the program 
//    that it is necessary to do some initialization work.  An input magnitude
//    of 2 lets the program know that that initialization can be skipped, 
//    and that useful information was computed earlier.)
//
//    The routine returns with all the information needed to continue
//    the integration.  If the integration reached TOUT, the user need only
//    define a new TOUT and call again.  In the one-step integrator
//    mode, returning with FLAG = -2, the user must keep in mind that 
//    each step taken is in the direction of the current TOUT.  Upon 
//    reaching TOUT, indicated by the output value of FLAG switching to 2,
//    the user must define a new TOUT and reset FLAG to -2 to continue 
//    in the one-step integrator mode.
//
//    In some cases, an error or difficulty occurs during a call.  In that case,
//    the output value of FLAG is used to indicate that there is a problem
//    that the user must address.  These values include:
//
//    * 3, integration was not completed because the input value of RELERR, the 
//      relative error tolerance, was too small.  RELERR has been increased 
//      appropriately for continuing.  If the user accepts the output value of
//      RELERR, then simply reset FLAG to 2 and continue.
//
//    * 4, integration was not completed because more than MAXNFE derivative 
//      evaluations were needed.  This is approximately (MAXNFE/6) steps.
//      The user may continue by simply calling again.  The function counter 
//      will be reset to 0, and another MAXNFE function evaluations are allowed.
//
//    * 5, integration was not completed because the solution vanished, 
//      making a pure relative error test impossible.  The user must use 
//      a non-zero ABSERR to continue.  Using the one-step integration mode 
//      for one step is a good way to proceed.
//
//    * 6, integration was not completed because the requested accuracy 
//      could not be achieved, even using the smallest allowable stepsize. 
//      The user must increase the error tolerances ABSERR or RELERR before
//      continuing.  It is also necessary to reset FLAG to 2 (or -2 when 
//      the one-step integration mode is being used).  The occurrence of 
//      FLAG = 6 indicates a trouble spot.  The solution is changing 
//      rapidly, or a singularity may be present.  It often is inadvisable 
//      to continue.
//
//    * 7, it is likely that this routine is inefficient for solving
//      this problem.  Too much output is restricting the natural stepsize
//      choice.  The user should use the one-step integration mode with 
//      the stepsize determined by the code.  If the user insists upon 
//      continuing the integration, reset FLAG to 2 before calling 
//      again.  Otherwise, execution will be terminated.
//
//    * 8, invalid input parameters, indicates one of the following:
//      NEQN <= 0;
//      T = TOUT and |FLAG| /= 1;
//      RELERR < 0 or ABSERR < 0;
//      FLAG == 0  or FLAG < -2 or 8 < FLAG.
//
//  Modified:
//
//    27 March 2004
//
//  Author:
//
//    H A Watts and L F Shampine,
//    Sandia Laboratories,
//    Albuquerque, New Mexico.
//
//  Reference:
//
//    E. Fehlberg,
//    Low-order Classical Runge-Kutta Formulas with Stepsize Control,
//    NASA Technical Report R-315.
//
//    L F Shampine, H A Watts, S Davenport,
//    Solving Non-stiff Ordinary Differential Equations - The State of the Art,
//    SIAM Review,
//    Volume 18, pages 376-411, 1976.
//
//  Parameters:
//
//    Input, external F, a user-supplied subroutine to evaluate the
//    derivatives Y'(T), of the form:
//
//      void f ( float t, float y[], float yp[] )
//
//    Input, int NEQN, the number of equations to be integrated.
//
//    Input/output, float Y[NEQN], the current solution vector at T.
//
//    Input/output, float YP[NEQN], the derivative of the current solution 
//    vector at T.  The user should not set or alter this information!
//
//    Input/output, float *T, the current value of the independent variable.
//
//    Input, float TOUT, the output point at which solution is desired.  
//    TOUT = T is allowed on the first call only, in which case the routine
//    returns with FLAG = 2 if continuation is possible.
//
//    Input, float *RELERR, ABSERR, the relative and absolute error tolerances
//    for the local error test.  At each step the code requires:
//      abs ( local error ) <= RELERR * abs ( Y ) + ABSERR
//    for each component of the local error and the solution vector Y.
//    RELERR cannot be "too small".  If the routine believes RELERR has been
//    set too small, it will reset RELERR to an acceptable value and return
//    immediately for user action.
//
//    Input, int FLAG, indicator for status of integration. On the first call, 
//    set FLAG to +1 for normal use, or to -1 for single step mode.  On 
//    subsequent continuation steps, FLAG should be +2, or -2 for single 
//    step mode.
//
//    Output, int RKF45_S, indicator for status of integration.  A value of 2 
//    or -2 indicates normal progress, while any other value indicates a 
//    problem that should be addressed.
//
{
# define MAXNFE 3000

  static float abserr_save = -1.0E+00;
  float ae;
  float dt;
  float ee;
  float eeoet;
  float eps;
  float esttol;
  float et;
  float *f1;
  float *f2;
  float *f3;
  float *f4;
  float *f5;
  static int flag_save = -1000;
  static float h = -1.0E+00;
  bool hfaild;
  float hmin;
  int i;
  static int init = -1000;
  int k;
  static int kflag = -1000;
  static int kop = -1;
  int mflag;
  static int nfe = -1;
  bool output;
  float relerr_min;
  static float relerr_save = -1.0E+00;
  static float remin = 1.0E-12;
  float s;
  float scale;
  float tol;
  float toln;
  float ypk;
//
//  Check the input parameters.
//
  eps = r_epsilon ( );

  if ( neqn < 1 )
  {
    return 8;
  }

  if ( (*relerr) < 0.0E+00 )
  {
    return 8;
  }

  if ( abserr < 0.0E+00 )
  {
    return 8;
  }

  if ( flag == 0 || 8 < flag  || flag < -2 )
  {
    return 8;
  }

  mflag = abs ( flag );
//
//  Is this a continuation call?
//
  if ( mflag != 1 )
  {
    if ( *t == tout && kflag != 3 )
    {
      return 8;
    }
//
//  FLAG = -2 or +2:
//
    if ( mflag == 2 )
    {
      if ( kflag == 3 )
      {
        flag = flag_save;
        mflag = abs ( flag );
      }
      else if ( init == 0 )
      {
        flag = flag_save;
      }
      else if ( kflag == 4 )
      {
        nfe = 0;
      }
      else if ( kflag == 5 && abserr == 0.0E+00 )
      {
        exit ( 1 );
      }
      else if ( kflag == 6 && (*relerr) <= relerr_save && abserr <= abserr_save )
      {
        exit ( 1 );
      }
    }
//
//  FLAG = 3, 4, 5, 6, 7 or 8.
//
    else
    {
      if ( flag == 3 )
      {
        flag = flag_save;
        if ( kflag == 3 )
        {
          mflag = abs ( flag );
        }
      }
      else if ( flag == 4 )
      {
        nfe = 0;
        flag = flag_save;
        if ( kflag == 3 )
        {
          mflag = abs ( flag );
        }
      }
      else if ( flag == 5 && 0.0E+00 < abserr )
      {
        flag = flag_save;
        if ( kflag == 3 )
        {
          mflag = abs ( flag );
        }
      }
//
//  Integration cannot be continued because the user did not respond to
//  the instructions pertaining to FLAG = 5, 6, 7 or 8.
//
      else
      {
        exit ( 1 );
      }
    }
  }
//
//  Save the input value of FLAG.  
//  Set the continuation flag KFLAG for subsequent input checking.
//
  flag_save = flag;
  kflag = 0;
//
//  Save RELERR and ABSERR for checking input on subsequent calls.
//
  relerr_save = (*relerr);
  abserr_save = abserr;
//
//  Restrict the relative error tolerance to be at least 
//
//    2*EPS+REMIN 
//
//  to avoid limiting precision difficulties arising from impossible 
//  accuracy requests.
//
  relerr_min = 2.0E+00 * r_epsilon ( ) + remin;
//
//  Is the relative error tolerance too small?
//
  if ( (*relerr) < relerr_min )
  {
    (*relerr) = relerr_min;
    kflag = 3;
    return 3;
  }

  dt = tout - *t;
//
//  Initialization:
//
//  Set the initialization completion indicator, INIT;
//  set the indicator for too many output points, KOP;
//  evaluate the initial derivatives
//  set the counter for function evaluations, NFE;
//  estimate the starting stepsize.
//
  f1 = new float[neqn];
  f2 = new float[neqn];
  f3 = new float[neqn];
  f4 = new float[neqn];
  f5 = new float[neqn];

  if ( mflag == 1 )
  {
    init = 0;
    kop = 0;
    f ( *t, y, yp );
    nfe = 1;

    if ( *t == tout )
    {
      return 2;
    }

  }

  if ( init == 0 )