RKF45method.cpp

//
//#####################################################################
//                          RKF45method.cpp 
//#####################################################################
//
// Class RKF45method implementation 
//
// For Math 157
//
//#####################################################################
// Chris Anderson (C) UCLA                             Jan. 22, 2009
//#####################################################################
//
#include "RKF45method.h"
#include "fortranODEfunction.h"

RKF45method::RKF45method() 
{
    errorTolerance = 1.0e-02;
    errorCode      = 0;
}

void RKF45method::setErrorTolerance(double tolerance)
{
    errorTolerance  = tolerance ;
}

double RKF45method::getErrorTolerance() const
{
    return errorTolerance;
}

int RKF45method::getErrorCode() const
{
    return errorCode;
}

//
// Prototype for the external ftc translation of rk45. 
//
extern "C" int rkf45_(int (*f)(double*, double*, double*), long *neqn, 
double *y, double*t, double *tout, double *relerr, 
double *abserr, long *iflag, double *work, long *iwork);

/**
@arg yk     : DoubleVector containing solution at time t = tk
@arg dt     : Timestep
@arg Fptr   : Pointer to to a class that extends the VectorFuntion class
\returns     
yk is updated with the solution value at t = tk + dt. The return value is 0 if advance is successful; 
otherwise the return value is the error code iflag from the rkf45 method. The value of 
iflag can also be obtained using the getErrorCode() member function. 
The following description is taken from the rkf45.f source code file. <p>

<pre>
       iflag = 3 -- integration was not completed because relative error
                    tolerance was too small. relerr has been increased
                    appropriately for continuing.
             = 4 -- integration was not completed because more than
                    3000 derivative evaluations were needed. this
                    is approximately 500 steps.
             = 5 -- integration was not completed because solution
                    vanished making a pure relative error test
                    impossible. must use 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 requested
                    accuracy could not be achieved using smallest
                    allowable stepsize. user must increase the error
                    tolerance before continued integration can be
                   attempted.
             = 7 -- it is likely that rkf45 is inefficient for solving
                    this problem. too much output is restricting the
                    natural stepsize choice. use the one-step integrator
                    mode.
             = 8 -- invalid input parameters
                    this indicator occurs if any of the following is
                    satisfied -   neqn .le. 0
                                  t=tout  and  iflag .ne. +1 or -1
                                  relerr or abserr .lt. 0.
                                  iflag .eq. 0  or  .lt. -2  or  .gt. 8
</pre>

The errorTolerance specified is used to set both the absolute and 
relative error tolerence values. As dicussed in the comments in
original rkf45.f, <a href="..\rkf45Comments.f">rkf45Comments.f</a>,
the advance() routine should not be used with relative error control 
smaller than about 1.e-8.
*/

int RKF45method::advance(DoubleVector& yk, double dt, VectorFunction* Fptr)
{
    long i;

    long dim = yk.getSize();
//
//  Use the fortranODEfunction class to provide 
//  a function pointer interface to the ODE specified by Fptr. 
//
    fortranODEfunction::setDimension(dim);
    fortranODEfunction::setFunction(Fptr);
    int (*f)(double*, double*, double*) = fortranODEfunction::f;
//
//  Allocate local temporary arrays
//
    double* y    = new double[dim];
    double* work = new double[4 + 6*dim];
    long*  iwork = new long[10];
//
//  Specify parameters and initial condition
//
    long neqn     = dim;
    double t      = 0.0;        // We assume an autonomous system; so always
    double tout   = dt;         // start at t = 0 and go to t = dt.
    double relerr = errorTolerance;
    double abserr = errorTolerance;
    long   iflag  = 1;

    for(i = 0; i < dim; i++)
    {
      y[i] =  yk(i);
    }
//
//  Call rkf45 method
//
    rkf45_(f,&neqn,y,&t,&tout,&relerr,&abserr,&iflag,work,iwork);
// 
//  Check for error in integration 
//
    if(iflag != 2)
    {
    // clean up 
    delete [] work;  delete [] iwork;  delete [] y;

    // Set error code to iflag and return iflag value 
    
    errorCode = iflag;
    return iflag;
    }
//
//  Call initialize so that ykp1 is guaranteed non-null.
// 
    
    for(i = 0; i < dim; i++)
    {
    yk(i) = y[i];
    }
    
    // clean up 

    delete [] work;  delete [] iwork; delete [] y;
    //
    // Set error code to 0 and return 0 for success 
    // Note: this code value is not the iflag value.
    //
    errorCode = 0;
    return 0;
}

---