NmOdeInt.cxx

Go to the documentation of this file.

// $Id: NmOdeInt.cxx,v 1.4 2005/03/03 21:48:29 musser Exp $
//
// A class to implement the integrations of ordinary differential
// equations.
//
// messier@huhepl.harvard.edu
#include "NumericalMethods/NmOdeInt.h"
#include <cassert>
#include <cmath>
#include "MessageService/MsgService.h"
#include "NumericalMethods/NmOdeDerivFunc.h"
#include "NumericalMethods/NmOdeStepperRKCK.h"
CVSID("$Id: NmOdeInt.cxx,v 1.4 2005/03/03 21:48:29 musser Exp $");

//......................................................................

NmOdeInt::NmOdeInt() :
  fXstart(0.0),
  fXend(0.0),
  fNvar(0),     // added by RWH
  fYstart(0),
  fYend(0),
  fYscale(0),
  fNmaxSteps(10000),
  fNok(0),
  fNbad(0),
  fEps(1.0e-6),
  fStepSize(0.0),
  fMinStepSize(0.0),
  fDerivs(0),
  fStepper(new NmOdeStepperRKCK), // Delete in d'ctor
  fDefaultStepper(true),
  fUseDefaultScales(true)
{ }

//......................................................................

NmOdeInt::NmOdeInt(const double* y1, 
                   int n, 
                   double x1, 
                   double x2,
                   double acc, 
                   double step, 
                   double minStep,
                   NmOdeDerivFunc& f) :
  fXstart(x1),
  fXend(x2),
  fNvar(n),
  fYstart(new double[fNvar]), // delete in d'ctor
  fYend(  new double[fNvar]), // delete in d'ctor
  fYscale(new double[fNvar]), // delete in d'ctor
  fNmaxSteps(10000),
  fNok(0),
  fNbad(0),
  fEps(acc),
  fStepSize(step),
  fMinStepSize(minStep),
  fDerivs(&f),
  fStepper(new NmOdeStepperRKCK), // delete in d'ctor
  fDefaultStepper(true),
  fUseDefaultScales(true)
{
//======================================================================
// Purpose: Build the class to perform the intergral
//
// Inputs: y1 - vector of initial dependent variable values
//          n - size of y1 vector, 
//         x1 - lower integration limit
//         x2 - upper integration limit
//        acc - fractional accurancy requested 
//       step - initial guess at a good step size for the interval
//    minStep - smallest step allowed on the interval
//          f - functor used to compute derivatives
//======================================================================
  // Initialize the start and end variable vectors
  for (int i=0; i<fNvar; ++i) {
    fYstart[i] = y1[i];
    fYend[i]   = 0.0;
  }
}

//......................................................................

NmOdeInt::~NmOdeInt() 
{
//======================================================================
// Purpose: Free memory allocated by class
//======================================================================
  delete []fYstart; fYstart = 0;
  delete []fYend;   fYend = 0;
  delete []fYscale; fYscale = 0;
  if (fDefaultStepper) { delete fStepper; fStepper = 0; }
}

//......................................................................

void NmOdeInt::SetMaxSteps(int n) 
{
//======================================================================
// Purpose: Set the maximum number of step to attempt in the interval
//======================================================================
  assert(n>1);
  fNmaxSteps = n;
}

//......................................................................

void NmOdeInt::SetNvar(int n) 
{
//======================================================================
// Purpose: Change the size of the vector of dependent variables
//======================================================================
  assert(n>0);
  // If the initial and final states have been allocated, delete them
  if (fYstart) { delete []fYstart; fYstart = 0; }
  if (fYend)   { delete []fYend;   fYend   = 0; }
  if (fYscale) { delete []fYscale; fYscale = 0; }
  fNvar = n;
  fYstart = new double[fNvar];
  fYend   = new double[fNvar];
  fYscale = new double[fNvar];
}

//......................................................................

void NmOdeInt::SetInitialCondition(const double* y) 
{
//======================================================================
// Purpose: Set the starting values for the dependent variables
//
// Inputs: y - the vector of variable values
//======================================================================
  assert(fNvar>0);
  for (int i=0; i<fNvar; ++i) fYstart[i] = y[i];
}

//......................................................................

void NmOdeInt::SetScales(const double* s) 
{
//======================================================================
// Purpose: Set the typical sizes of the dependent variables
//
// Inputs: s - vector of typical sizes for the dependent variables
//======================================================================
  assert(fNvar>0);
  for (int i=0; i<fNvar; ++i) fYscale[i] = s[i];
  fUseDefaultScales = false;
}

//......................................................................

void NmOdeInt::SetStepper(NmOdeStepper& s) 
{
//======================================================================
// Purpose: Change the stepper class to use to perform the integral
//
// Inputs: s - the new stepper to use
//======================================================================
  if (fStepper && fDefaultStepper) { 
    delete fStepper; 
    fStepper = 0; 
  }
  fDefaultStepper = false; // No longer the owner of the stepper
  fStepper = &s;
}

//......................................................................

bool NmOdeInt::Integrate(double* yFinal) 
{
//======================================================================
// Purpose: Perform the integral
//
// Output: yFinal - the final values of the dependent variables
//
// Returns: true  - success
//          false - failure
//======================================================================
  register int nstp; // Loop variable
  register int i;    // Loop variable
//unused:  double xSave;
  double x;
  double hNext;
  double hDid;
  double h;
  double*    y = new double[fNvar]; 
  double* dydx = new double[fNvar];
  static const double tiny = 1.0e-30;

  // Check and set initial conditions and provide some defaults
  fNok = fNbad = 0;
  if (fEps == 0.0) fEps = 1.0e-4;
  if (fStepSize == 0.0) 
    fStepSize = sqrt(fEps)*fabs(fXstart-fXend);
  
  x = fXstart;
  if (fXend>fXstart) h =  fStepSize;
  else               h = -fStepSize;

  for (i=0; i<fNvar; ++i) y[i] = fYstart[i];
  for (nstp = 0; nstp<fNmaxSteps; ++nstp) {
    (*fDerivs)(x,y,dydx);
    
    // Set the scale for the variables
    if (fUseDefaultScales) {
      for (i=0; i<fNvar; ++i) {
        fYscale[i] = fabs(y[i])+fabs(dydx[i]*h)+tiny;
      }
    }
    
    // Check if step size over shoots end point
    if ((x+h-fXend)*(x+h-fXstart) > 0.0) h = fXend-x;
    // Take the step
    fStepper->Step(y, dydx, fNvar, &x, h, fEps, fYscale, 
                   &hDid, &hNext, *fDerivs);
    if (hDid == h) fNok++;
    else           fNbad++;
    
    // Check if we're done
    if (this->Done(x,y)) {
      // Finish up, free memory, and return
      for (i=0; i<fNvar; ++i) {
        yFinal[i] = fYend[i] = y[i];
      }
      delete []y;    y = 0;    // deallocate temporary arrays
      delete []dydx; dydx = 0; // deallocate temporary arrays

      return true;
    }
    if (fabs(hNext)<=fMinStepSize) {
      MSG("Nm",Msg::kWarning) << 
        "Step size too small (" << hNext<< "<" << fMinStepSize << "\n";
      h = fMinStepSize;
    }
    else {
      h = hNext;
    }
  }
  // Only get here when we've exceeded the max. number of steps
  for (i=0; i<fNvar; ++i) {
    yFinal[i] = fYend[i] = y[i];
  }
  delete []y;    y = 0;    // deallocate temporary arrays
  delete []dydx; dydx = 0; // deallocate temporary arrays
  MSG("Nm",Msg::kError) 
    << "Maximum number of steps ("<< fNmaxSteps <<") exceeded\n";
  return false;
}

//......................................................................

bool NmOdeInt::Done(double x, const double* y) 
{
//======================================================================
// Purpose: Test if the integral is complete
//
// Inputs: x - value of integration variable
//       y - values of dependent variables
//
// Returns: true  - done
//          false - not done
//======================================================================
  return ((x-fXend)*(fXend-fXstart)>=0.0);
}

---