---

QualityControlCashKarpMover.cpp

Go to the documentation of this file.

// QualityControlCashKarpMover.cpp: implementation of the QualityControlCashKarpMover class.
//

#include "QualityControlCashKarpMover.h"

// Parameters for 5th order step

const double CQualityControlCashKarpMover::m_scdA[6] = {
        0.0, 0.2, 0.3, 0.6, 1.0, 0.875  };

const double CQualityControlCashKarpMover::m_scdB[6][5] = {
                { 0.0,                 0.0,                  0.0,                        0.0,                           0.0                     },
                { 1.0/5.0,                 0.0,              0.0,                        0.0,                           0.0                     },
                { 3.0/40.0,                9.0/40.0,     0.0,                    0.0,                           0.0                     },
                { 3.0/10.0,               -9.0/10.0,     6.0/5.0,            0.0,                               0.0                     },
                {-11.0/54.0,       5.0/2.0,         -70.0/27.0,      35.0/27.0,                 0.0                     },
                { 1631.0/55296.0,  175.0/512.0,  575.0/13824.0,  44275.0/110592.0,  253.0/4096.0}  };

const double CQualityControlCashKarpMover::m_scdC[6] = {
        37.0/378.0, 0.0, 250.0/621.0, 125.0/594.0, 0.0, 512.0/1771.0 };

const double CQualityControlCashKarpMover::m_scdCStar[6] = {            
        2825.0/27648.0, 0.0, 18575.0/48384.0, 13525.0/55296.0, 277.0/14336.0, 0.25  };

const double CQualityControlCashKarpMover::m_scdERRCON = 1.89e-04;

const double CQualityControlCashKarpMover::m_scdPGROW = -0.2;

const double CQualityControlCashKarpMover::m_scdPSHRINK = -0.25;

const double CQualityControlCashKarpMover::m_scdSAFETY = 0.9;

const double CQualityControlCashKarpMover::m_scdTINY = 1.0e-30;

// Construction/Destruction

CQualityControlCashKarpMover::CQualityControlCashKarpMover(double frequency, double stepSize, double eps)
:CDifferentialEquationMover(frequency,stepSize)
{
        m_dLastStep = stepSize;
        m_dNextStep = stepSize;
        m_dInitialStepSize = stepSize;
        m_dEps = eps;
}

CQualityControlCashKarpMover::~CQualityControlCashKarpMover()

{
        
}


void CQualityControlCashKarpMover::Move(double xInitial, double xFinal, ReactionNetwork *pReactionNetwork)
{
        int i;
        // do a couple of checks on the integration interval
        double xInterval = xFinal-xInitial;
        // set the member pointer equal to the pointer passed in
        m_pReactionNetwork = pReactionNetwork;
        // if there's no moving to be done and xI = xF = integer, take data there
        if(this->MoveTimeIsZero(xInterval)) 
        {
                if(xInitial > (int) xInitial)
                {
                }
                else
                {
                        m_iCount = (int) xInitial;
                        m_iCount++;
                        Notify();
                }
                // take a data point 
                return;
        }
        if(m_dInitialStepSize > xInterval) {m_dInitialStepSize = xInterval;}

        m_dTime = xInitial;
        double nextTime;
        // set up the count variable
        m_iCount = (int) xInitial;
        m_iCount++;
        nextTime = m_iCount;

        m_dStepSize = m_dInitialStepSize;
        m_dLastStep = m_dInitialStepSize;
        m_dNextStep = m_dInitialStepSize;
        
        int nChem = m_pReactionNetwork->GetNumberOfChemicals();
        double *chem = new double[nChem];
        double *dChemdt = new double[nChem];
        m_pdScale = new double[nChem];
        m_pdError = new double[nChem];

        // fill chem[] with starting chemical concentrations
        for(i = 0; i < nChem; i++)
        {
                chem[i] = m_pReactionNetwork->GetChemical(i)->GetAmount();
        }

        while(m_dTime < xFinal)
        {
                // run until xFinal or nextTime (nextTime might be > xFinal)
                while( m_dTime < __min(xFinal,nextTime) )
                {
                        // get initial derivatives
                        ComputeDerivatives(chem,dChemdt);
                        // compute scaling vector
                        // __max is to ensure chemicals that start at ~ 0 don't get
                        // huge scaling factors, which screws up the stepper
                        for(i = 0; i < nChem; i++)
                        {
                                m_pdScale[i] = fabs(chem[i]) + fabs(m_dStepSize*dChemdt[i]) + m_scdTINY;
                        }
                        
                        // don't run over - nextTime might be > xFinal
                        if( (((m_dTime + m_dNextStep) - nextTime) > 0.0) || (((m_dTime + m_dNextStep) - xFinal) > 0.0) )        
                        {
                                m_dNextStep = __min(nextTime-m_dTime,xFinal-m_dTime);
                        }

                        // call the stepper
                        Stepper(nChem,chem,dChemdt);

                        
                        // ensure synchrony with chem and network
                        for(int i = 0; i < nChem; i++)
                        {
                                m_pReactionNetwork->GetChemical(i)->SetAmount(chem[i]);
                        }
                        
                        // record largest/smallest/numsteps info
                        if(m_dLastStep < this->GetMinStepSize()) {this->SetMinStepSize(m_dLastStep);}
                        if(m_dLastStep > this->GetMaxStepSize()) {this->SetMaxStepSize(m_dLastStep);}
                        this->IncrementTotalSteps();
                }
                
                if(m_dTime == nextTime)
                {
                        m_iCount++;
                        Notify();
                        nextTime += 1.0;
                }
        }

        delete [] chem;
        delete [] dChemdt;
        delete [] m_pdScale;
        delete [] m_pdError;
}

// The stepper adjusts the stepsize according to the error and takes
// a single 5th order Cash-Karp step

void CQualityControlCashKarpMover::Stepper(int nRHS, double *y, double *dydt)
{
        int i;
        double *yTemp = new double[nRHS];
        double maxError;
        
        m_dStepSize = m_dNextStep;

        while(1)
        {
                for(i = 0; i < nRHS; i++)
                {
                        yTemp[i] = y[i];
                }
                
                // try a step
                CashKarpStep(nRHS,yTemp,dydt);
                // evaluate accuracy
                maxError = 0.0;
                for(i = 0; i < nRHS; i++)
                {
                        maxError = __max(maxError,fabs(m_pdError[i]/m_pdScale[i]));
                }
                maxError = maxError/m_dEps;

                // step successful; break out
                if(maxError <= 1.0) {break;}
                // step unsuccessful, make stepsize smaller
                double tempStep = m_scdSAFETY*m_dStepSize*pow(maxError,m_scdPSHRINK);
                
                if(m_dStepSize >= 0.0)
                {
                        m_dStepSize = __max(tempStep,0.1*m_dStepSize);
                }
                else
                {
                        m_dStepSize = __min(tempStep,0.1*m_dStepSize);
                }
                if((m_dTime + m_dStepSize) == m_dTime)
                {
                        // step extremely tiny, should throw an exception
                        cout << "Stepsize too small in QualityControlCashKarpMover::Stepper()" << endl;
                        exit(1);
                }
        }
        // we accepted the step, so m_dLastStep gets updated
        m_dLastStep = m_dStepSize;
        if(maxError > m_scdERRCON) 
        {
                m_dNextStep = m_scdSAFETY*m_dStepSize*pow(maxError,m_scdPGROW);
        }
        else 
        {
                m_dNextStep = 5.0*m_dStepSize;
        }
        // increment m_dTime
        m_dTime += m_dLastStep;
        // accept the solution
        for(int i = 0; i < nRHS; i++)
        {
                y[i] = yTemp[i];
        }

        delete [] yTemp;
}

// Takes one 5th order Runge-Kutta-Fehlberg step using parameters 
// from Cash and Karp.  The new solution is returned in y[].  The 
// error in the solution is estimated as the difference between 4th
// and 5th order methods.

void CQualityControlCashKarpMover::CashKarpStep(int nRHS, double *y, double *dydt)
{
        int i;
        double *dy2 = new double[nRHS];
        double *dy3 = new double[nRHS];
        double *dy4 = new double[nRHS];
        double *dy5 = new double[nRHS];
        double *dy6 = new double[nRHS];
        double *yTemp = new double[nRHS];

        // first step
        for(i = 0; i < nRHS; i++)
        {
                yTemp[i] = y[i] + m_scdB[1][0]*m_dStepSize*dydt[i];
        }
        ComputeDerivatives(yTemp,dy2);
        
        // second step
        for(i = 0; i < nRHS; i++)
        {
                yTemp[i] = y[i] + m_dStepSize*(m_scdB[2][0]*dydt[i] + m_scdB[2][1]*dy2[i]);
        }
        ComputeDerivatives(yTemp,dy3);

        // third step
        for(i = 0; i < nRHS; i++)
        {
                yTemp[i] = y[i] + m_dStepSize*(m_scdB[3][0]*dydt[i] + m_scdB[3][1]*dy2[i] + m_scdB[3][2]*dy3[i]);
        }
        ComputeDerivatives(yTemp,dy4);

        // fourth step
        for(i = 0; i < nRHS; i++)
        {
                yTemp[i] = y[i] + m_dStepSize*(m_scdB[4][0]*dydt[i] + m_scdB[4][1]*dy2[i] + m_scdB[4][2]*dy3[i] + m_scdB[4][3]*dy4[i]);
        }
        ComputeDerivatives(yTemp,dy5);

        // fifth step
        for(i = 0; i < nRHS; i++)
        {
                yTemp[i] = y[i] + m_dStepSize*(m_scdB[5][0]*dydt[i] + m_scdB[5][1]*dy2[i] + m_scdB[5][2]*dy3[i] + m_scdB[5][3]*dy4[i] + m_scdB[5][4]*dy5[i]);
        }
        ComputeDerivatives(yTemp,dy6);

        // form new solution and estimate error
        // error estimate is the difference between 4th and 5th order methods
        for(i = 0; i < nRHS; i++)
        {
                y[i] = y[i] + m_dStepSize*(m_scdC[0]*dydt[i] + m_scdC[2]*dy3[i] + m_scdC[3]*dy4[i] + m_scdC[5]*dy6[i]);
                // this form of the calculation is to try to minimize roundoff error
                double temp1 = m_dStepSize*(m_scdC[0] - m_scdCStar[0])*dydt[i];
                double temp2 = m_dStepSize*(m_scdC[2] - m_scdCStar[2])*dy3[i];
                double temp3 = m_dStepSize*(m_scdC[3] - m_scdCStar[3])*dy4[i];
                double temp4 = m_dStepSize*(m_scdC[4] - m_scdCStar[4])*dy5[i];
                double temp5 = m_dStepSize*(m_scdC[5] - m_scdCStar[5])*dy6[i];
                m_pdError[i] = temp1 + temp2 + temp3 + temp4 + temp5;
        }
        

        delete [] dy2;
        delete [] dy3;
        delete [] dy4;
        delete [] dy5;
        delete [] dy6;
        delete [] yTemp;
}

---