---

RK4Mover.cpp

Go to the documentation of this file.

// RK4Mover.cpp: implementation of the CRK4Mover class.
//

#include "RK4Mover.h"

// Construction/Destruction

CRK4Mover::CRK4Mover()
{

}

CRK4Mover::CRK4Mover(double frequency, double stepSize)
:CDifferentialEquationMover(frequency,stepSize)
{

}

CRK4Mover::~CRK4Mover()
{

}

void CRK4Mover::Move(double xInitial, double xFinal, ReactionNetwork *pReactionNetwork)
{
        // do a couple of checks on the integration interval
        double xInterval = xFinal-xInitial;
        if(this->MoveTimeIsZero(xInterval)) {return;}
        if(m_dStepSize > xInterval) {m_dStepSize = xInterval;}

        // set the member pointer to the network pointer passed in
        this->m_pReactionNetwork = pReactionNetwork;

        int dataCnt=0;
        int nChem = m_pReactionNetwork->GetNumberOfChemicals();
        double *chem = new double[nChem];
        double *dChemdt = new double[nChem];
        int freq = 1.0/m_dStepSize;
        
        m_dTime = xInitial;
        // may be unsafe if xInitial has a fractional part! Most likely the storage of time 
        // series data will be altered before this is a problem.
        m_iCount = (int) xInitial;
        double nextTime = m_dTime;

        for (int cNum = 0; cNum < nChem; cNum++)
        {
                chem[cNum] = m_pReactionNetwork->GetChemical(cNum)->GetAmount();
        }       

        while(m_dTime < xFinal)
        {
                while(m_dTime < nextTime)
                {
                        ComputeDerivatives(chem,dChemdt);
                        FourthOrderStep(nChem,chem,dChemdt);

                        m_dTime += m_dStepSize;
                }

                // before notification, make sure chemical conc. are current
                for(int i = 0; i < nChem; i++)
                {
                        m_pReactionNetwork->GetChemical(i)->SetAmount(chem[i]);
                }

                m_iCount++;
                Notify();
                nextTime += 1.0;
        }

        // free memory

        delete [] chem;
        delete [] dChemdt;
}

// Takes a single Runge-Kutta 4th order step of size stepSize (whatever
// its current value).  Upon exit, the new values of the variables are
// stored in y[]

void CRK4Mover::FourthOrderStep(int nRHS, double *y, double *dydt)
{
        int i;
        double *yTemp = new double[nRHS];
        double *dydtTemp1 = new double[nRHS];
        double *dydtTemp2 = new double[nRHS];
        
        double halfStep = 0.5*m_dStepSize;
        double stepOverSix = m_dStepSize/6.0;
        
        for(i = 0; i < nRHS; i++)
        {
                yTemp[i] = y[i] + halfStep*dydt[i];
        }
        ComputeDerivatives(yTemp,dydtTemp1);

        for(i = 0; i < nRHS; i++)
        {
                yTemp[i] = y[i] + halfStep*dydtTemp1[i];
        }
        ComputeDerivatives(yTemp,dydtTemp2);

        for(i = 0; i < nRHS; i++)
        {
                yTemp[i] = y[i] + m_dStepSize*dydtTemp2[i];
                dydtTemp2[i] += dydtTemp1[i];
        }
        ComputeDerivatives(yTemp,dydtTemp1);
        
        for(i = 0; i < nRHS; i++)
        {
                y[i] = y[i] + stepOverSix*(dydt[i] + dydtTemp1[i] + 2.0*dydtTemp2[i]);
        }

        // free memory
        delete [] yTemp;
        delete [] dydtTemp1;
        delete [] dydtTemp2;
}

---