rkf45.c

Runge-Kutta-Fehlberg method

# include <cstdlib>
# include <iostream>
# include <iomanip>
# include <cmath>
# include <ctime>

using namespace std;

# include "rkf45.H"

//******************************************************************************

double d_epsilon ( void )

//******************************************************************************
//
//  Purpose:
//
//    D_EPSILON returns the round off unit for double precision arithmetic.
//
//  Discussion:
//
//    D_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 D_EPSILON, the floating point round-off unit.
//
{
  float r;

  r = 1.0E+00;

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

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

double d_max ( double x, double y )

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

double d_min ( double x, double y )

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

double d_sign ( double x )

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

void fehl_d ( void f ( double t, double y[], double yp[] ), int neqn, 
  double y[], double t, double h, double yp[], double f1[], double f2[], double f3[], 
  double f4[], double f5[], double s[] )

//******************************************************************************
//
//  Purpose:
//
//    FEHL_D takes one Fehlberg fourth-fifth order step (double precision).
//
//  Discussion:
//
//    This routine integrates a system of NEQN first order ordinary differential
//    equations of the form
//      dY(i)/dT = F(T,Y(1:NEQN))
//    where the initial values Y and the initial derivatives
//    YP are specified at the starting point T.
//
//    The routine advances the solution over the fixed step H and returns
//    the fifth order (sixth order accurate locally) solution
//    approximation at T+H in array S.
//
//    The formulas have been grouped to control loss of significance.
//    The routine should be called with an H not smaller than 13 units of
//    roundoff in T so that the various independent arguments can be
//    distinguished.
//
//  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 ( double t, double y[], double yp[] )
//
//    Input, int NEQN, the number of equations to be integrated.
//
//    Input, double Y[NEQN], the current value of the dependent variable.
//
//    Input, double T, the current value of the independent variable.
//
//    Input, double H, the step size to take.
//
//    Input, double YP[NEQN], the current value of the derivative of the
//    dependent variable.
//
//    Output, double F1[NEQN], F2[NEQN], F3[NEQN], F4[NEQN], F5[NEQN], derivative
//    values needed for the computation.
//
//    Output, double S[NEQN], the estimate of the solution at T+H.
//
{
  double ch;
  int i;

  ch = h / 4.0E+00;

  for ( i = 0; i < neqn; i++ )
  {
    f5[i] = y[i] + ch * yp[i];
  }

  f ( t + ch, f5, f1 );

  ch = 3.0E+00 * h / 32.0E+00;

  for ( i = 0; i < neqn; i++ )
  {
    f5[i] = y[i] + ch * ( yp[i] + 3.0E+00 * f1[i] );
  }

  f ( t + 3.0E+00 * h / 8.0E+00, f5, f2 );

  ch = h / 2197.0E+00;

  for ( i = 0; i < neqn; i++ )
  {
    f5[i] = y[i] + ch * 
    ( 1932.0E+00 * yp[i] 
    + ( 7296.0E+00 * f2[i] - 7200.0E+00 * f1[i] ) 
    );
  }

  f ( t + 12.0E+00 * h / 13.0E+00, f5, f3 );

  ch = h / 4104.0E+00;

  for ( i = 0; i < neqn; i++ )
  {
    f5[i] = y[i] + ch * 
    ( 
      ( 8341.0E+00 * yp[i] - 845.0E+00 * f3[i] ) 
    + ( 29440.0E+00 * f2[i] - 32832.0E+00 * f1[i] ) 
    );
  }

  f ( t + h, f5, f4 );

  ch = h / 20520.0E+00;

  for ( i = 0; i < neqn; i++ )
  {
    f1[i] = y[i] + ch * 
    ( 
      ( -6080.0E+00 * yp[i] 
      + ( 9295.0E+00 * f3[i] - 5643.0E+00 * f4[i] ) 
      ) 
    + ( 41040.0E+00 * f1[i] - 28352.0E+00 * f2[i] ) 
    );
  }

  f ( t + h / 2.0E+00, f1, f5 );
//
//  Ready to compute the approximate solution at T+H.
//
  ch = h / 7618050.0E+00;

  for ( i = 0; i < neqn; i++ )
  {
    s[i] = y[i] + ch * 
    ( 
      ( 902880.0E+00 * yp[i] 
      + ( 3855735.0E+00 * f3[i] - 1371249.0E+00 * f4[i] ) ) 
    + ( 3953664.0E+00 * f2[i] + 277020.0E+00 * f5[i] ) 
    );
  }

  return;
}
//******************************************************************************

void fehl_s ( void f ( float t, float y[], float yp[] ), int neqn, 
  float y[], float t, float h, float yp[], float f1[], float f2[], float f3[], 
  float f4[], float f5[], float s[] )

//******************************************************************************
//
//  Purpose:
//
//    FEHL_S takes one Fehlberg fourth-fifth order step (single precision).
//
//  Discussion:
//
//    This routine integrates a system of NEQN first order ordinary differential
//    equations of the form
//      dY(i)/dT = F(T,Y(1:NEQN))
//    where the initial values Y and the initial derivatives
//    YP are specified at the starting point T.
//
//    The routine advances the solution over the fixed step H and returns
//    the fifth order (sixth order accurate locally) solution
//    approximation at T+H in array S.
//
//    The formulas have been grouped to control loss of significance.
//    The routine should be called with an H not smaller than 13 units of
//    roundoff in T so that the various independent arguments can be
//    distinguished.
//
//  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 ( double t, double y[], double yp[] )
//
//    Input, int NEQN, the number of equations to be integrated.
//
//    Input, float Y[NEQN], the current value of the dependent variable.
//
//    Input, float T, the current value of the independent variable.
//
//    Input, float H, the step size to take.
//
//    Input, float YP[NEQN], the current value of the derivative of the
//    dependent variable.
//
//    Output, float F1[NEQN], F2[NEQN], F3[NEQN], F4[NEQN], F5[NEQN], derivative
//    values needed for the computation.
//
//    Output, float S[NEQN], the estimate of the solution at T+H.
//
{
  float ch;
  int i;

  ch = h / 4.0E+00;

  for ( i = 0; i < neqn; i++ )
  {
    f5[i] = y[i] + ch * yp[i];
  }

  f ( t + ch, f5, f1 );

  ch = 3.0E+00 * h / 32.0E+00;

  for ( i = 0; i < neqn; i++ )
  {
    f5[i] = y[i] + ch * ( yp[i] + 3.0E+00 * f1[i] );
  }

  f ( t + 3.0E+00 * h / 8.0E+00, f5, f2 );

  ch = h / 2197.0E+00;

  for ( i = 0; i < neqn; i++ )
  {
    f5[i] = y[i] + ch * 
    ( 1932.0E+00 * yp[i] 
    + ( 7296.0E+00 * f2[i] - 7200.0E+00 * f1[i] ) 
    );
  }

  f ( t + 12.0E+00 * h / 13.0E+00, f5, f3 );

  ch = h / 4104.0E+00;

  for ( i = 0; i < neqn; i++ )
  {
    f5[i] = y[i] + ch * 
    ( 
      ( 8341.0E+00 * yp[i] - 845.0E+00 * f3[i] ) 
    + ( 29440.0E+00 * f2[i] - 32832.0E+00 * f1[i] ) 
    );
  }

  f ( t + h, f5, f4 );

  ch = h / 20520.0E+00;

  for ( i = 0; i < neqn; i++ )
  {
    f1[i] = y[i] + ch * 
    ( 
      ( -6080.0E+00 * yp[i] 
      + ( 9295.0E+00 * f3[i] - 5643.0E+00 * f4[i] ) 
      ) 
    + ( 41040.0E+00 * f1[i] - 28352.0E+00 * f2[i] ) 
    );
  }

  f ( t + h / 2.0E+00, f1, f5 );
//
//  Ready to compute the approximate solution at T+H.
//
  ch = h / 7618050.0E+00;

  for ( i = 0; i < neqn; i++ )
  {
    s[i] = y[i] + ch * 
    ( 
      ( 902880.0E+00 * yp[i] 
      + ( 3855735.0E+00 * f3[i] - 1371249.0E+00 * f4[i] ) ) 
    + ( 3953664.0E+00 * f2[i] + 277020.0E+00 * f5[i] ) 
    );
  }

  return;
}
//******************************************************************************

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

//******************************************************************************
//
//  Purpose:
//
//    RKF45_D carries out the Runge-Kutta-Fehlberg method (double 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: