rkf45.c

Runge-Kutta-Fehlberg method

  {
    init = 1;
    h = fabs ( dt );
    toln = 0.0E+00;

    for ( k = 0; k < neqn; k++ )
    {
      tol = (*relerr) * fabs ( y[k] ) + abserr;
      if ( 0.0E+00 < tol )
      {
        toln = tol;
        ypk = fabs ( yp[k] );
        if ( tol < ypk * pow ( h, 5 ) )
        {
          h = ( float ) pow ( ( double ) ( tol / ypk ), 0.2E+00 );
        }
      }
    }

    if ( toln <= 0.0E+00 )
    {
      h = 0.0E+00;
    }

    h = r_max ( h, 26.0E+00 * eps * r_max ( fabs ( *t ), fabs ( dt ) ) );

    if ( flag < 0 )
    {
      flag_save = -2;
    }
    else
    {
      flag_save = 2;
    }
  }
//
//  Set stepsize for integration in the direction from T to TOUT.
//
  h =  r_sign ( dt ) * fabs ( h );
//
//  Test to see if too may output points are being requested.
//
  if ( 2.0E+00 * fabs ( dt ) <= fabs ( h ) )
  {
    kop = kop + 1;
  }
//
//  Unnecessary frequency of output.
//
  if ( kop == 100 )
  {
    kop = 0;
    delete [] f1;
    delete [] f2;
    delete [] f3;
    delete [] f4;
    delete [] f5;
    return 7;
  }
//
//  If we are too close to the output point, then simply extrapolate and return.
//
  if ( fabs ( dt ) <= 26.0E+00 * eps * fabs ( *t ) )
  {
    *t = tout;
    for ( i = 0; i < neqn; i++ )
    {
      y[i] = y[i] + dt * yp[i];
    }
    f ( *t, y, yp );
    nfe = nfe + 1;

    delete [] f1;
    delete [] f2;
    delete [] f3;
    delete [] f4;
    delete [] f5;
    return 2;
  }
//
//  Initialize the output point indicator.
//
  output = false;
//
//  To avoid premature underflow in the error tolerance function,
//  scale the error tolerances.
//
  scale = 2.0E+00 / (*relerr);
  ae = scale * abserr;
//
//  Step by step integration.
//
  for ( ; ; )
  {
    hfaild = false;
//
//  Set the smallest allowable stepsize.
//
    hmin = 26.0E+00 * eps * fabs ( *t );
//
//  Adjust the stepsize if necessary to hit the output point.
//
//  Look ahead two steps to avoid drastic changes in the stepsize and
//  thus lessen the impact of output points on the code.
//
    dt = tout - *t;

    if ( 2.0E+00 * fabs ( h ) <= fabs ( dt ) )
    {
    }
    else
//
//  Will the next successful step complete the integration to the output point?
//
    {
      if ( fabs ( dt ) <= fabs ( h ) )
      {
        output = true;
        h = dt;
      }
      else
      {
        h = 0.5E+00 * dt;
      }

    }
//
//  Here begins the core integrator for taking a single step.
//
//  The tolerances have been scaled to avoid premature underflow in
//  computing the error tolerance function ET.
//  To avoid problems with zero crossings, relative error is measured
//  using the average of the magnitudes of the solution at the
//  beginning and end of a step.
//  The error estimate formula has been grouped to control loss of
//  significance.
//
//  To distinguish the various arguments, H is not permitted
//  to become smaller than 26 units of roundoff in T.
//  Practical limits on the change in the stepsize are enforced to
//  smooth the stepsize selection process and to avoid excessive
//  chattering on problems having discontinuities.
//  To prevent unnecessary failures, the code uses 9/10 the stepsize
//  it estimates will succeed.
//
//  After a step failure, the stepsize is not allowed to increase for
//  the next attempted step.  This makes the code more efficient on
//  problems having discontinuities and more effective in general
//  since local extrapolation is being used and extra caution seems
//  warranted.
//
//  Test the number of derivative function evaluations.
//  If okay, try to advance the integration from T to T+H.
//
    for ( ; ; )
    {
//
//  Have we done too much work?
//
      if ( MAXNFE < nfe )
      {
        kflag = 4;
        delete [] f1;
        delete [] f2;
        delete [] f3;
        delete [] f4;
        delete [] f5;
        return 4;
      }
//
//  Advance an approximate solution over one step of length H.
//
      fehl_s ( 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 = r_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 / ( float ) pow ( ( double ) 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 / ( float ) pow ( ( double ) esttol, 0.2E+00 );
    }
    else
    {
      s = 5.0E+00;
    }

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

    h = r_sign ( h ) * r_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
}
//******************************************************************************

float r_epsilon ( void )

//******************************************************************************
//
//  Purpose:
//
//    R_EPSILON returns the round off unit for floating arithmetic.
//
//  Discussion:
//
//    R_EPSILON is a number R which is a power of 2 with the property that,
//    to the precision of the computer's arithmetic,
//      1 < 1 + R
//    but 
//      1 = ( 1 + R / 2 )
//
//  Modified:
//
//    06 May 2003
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, float R_EPSILON, the floating point round-off unit.
//
{
  float r;

  r = 1.0E+00;

  while ( 1.0E+00 < ( float ) ( 1.0E+00 + r )  )
  {
    r = r / 2.0E+00;
  }

  return ( 2.0E+00 * r );
}
//*********************************************************************

float r_max ( float x, float y )

//*********************************************************************
//
//  Purpose:
//
//    R_MAX returns the maximum of two real values.
//
//  Modified:
//
//    10 January 2002
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, float X, Y, the quantities to compare.
//
//    Output, float R_MAX, the maximum of X and Y.
//
{
  if ( y < x )
  {
    return x;
  } 
  else
  {
    return y;
  }
}
//*********************************************************************

float r_min ( float x, float y )

//*********************************************************************
//
//  Purpose:
//
//    R_MIN returns the minimum of two real values.
//
//  Modified:
//
//    09 May 2003
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, float X, Y, the quantities to compare.
//
//    Output, float R_MIN, the minimum of X and Y.
//
{
  if ( y < x )
  {
    return y;
  } 
  else
  {
    return x;
  }
}
//*********************************************************************

float r_sign ( float x )

//*********************************************************************
//
//  Purpose:
//
//    R_SIGN returns the sign of a real number.
//
//  Modified:
//
//    27 March 2004
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, float X, the number whose sign is desired.
//
//    Output, float R_SIGN, the sign of X.
//
{
  if ( x < 0.0E+00 )
  {
    return ( -1.0E+00 );
  } 
  else
  {
    return ( +1.0E+00 );
  }
}
//**********************************************************************

void timestamp ( void )

//**********************************************************************
//
//  Purpose:
//
//    TIMESTAMP prints the current YMDHMS date as a time stamp.
//
//  Example:
//
//    May 31 2001 09:45:54 AM
//
//  Modified:
//
//    24 September 2003
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    None
//
{
#define TIME_SIZE 40

  static char time_buffer[TIME_SIZE];
  const struct tm *tm;
  size_t len;
  time_t now;

  now = time ( NULL );
  tm = localtime ( &now );

  len = strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm );

  cout << time_buffer << "\n";

  return;
#undef TIME_SIZE
}