RungeKutta.h

Go to the documentation of this file.

#ifndef RUNGEKUTTA_H
#define RUNGEKUTTA_H

#include <Matrix.h>

typedef Vector(*ptr2RKFunc)(const double &_time, const Vector &_states, const Matrix &_params);

inline Matrix RungeKuttaSolve(Vector (*odeFunc)(const double &timeInput, const Vector &stateInput, const Matrix &constantsInput), const Vector &_initialConditions, const Matrix &_constants, const double &_timeInitial, const double &_timeFinal, const int &_numSteps)
{
        double h = (_timeFinal - _timeInitial) / _numSteps;
        double t = _timeInitial;
        int numEqs = _initialConditions.getIndexCount();
        Vector inputs(numEqs); inputs = _initialConditions;
        Vector K1(numEqs);
        Vector K2(numEqs);
        Vector K3(numEqs);
        Vector K4(numEqs);

        Matrix output(_numSteps + 1, numEqs + 1);

        output(MatrixIndexBase,MatrixIndexBase) = _timeInitial;
        output(MatrixIndexBase,_(MatrixIndexBase+1,MatrixIndexBase+numEqs)) = ~_initialConditions;

        Vector temp(numEqs);
        for (int ii = 1; ii <= _numSteps; ++ii)
        {
            K1 = h * odeFunc(t, inputs, _constants);
            temp = inputs + K1 / 2.0;
            K2 = h * odeFunc(t + h/2, temp, _constants);
            temp = inputs + K2 / 2.0;
            K3 = h * odeFunc(t + h/2, temp, _constants);
            temp = inputs + K3;
            K4 = h * odeFunc(t + h, temp, _constants);
            for (int jj = MatrixIndexBase; jj < MatrixIndexBase + numEqs; ++jj)
            {
                inputs(jj) += (K1(jj)
                                 + 2.0 * K2(jj)
                                 + 2.0 * K3(jj)
                                 + K4(jj)) / 6.0;                       
            }
            t = _timeInitial + ii * h;
            output(MatrixIndexBase + ii,MatrixIndexBase) = t;
            output(MatrixIndexBase + ii,_(MatrixIndexBase+1,MatrixIndexBase+numEqs)) = ~inputs; 
        }
        return output;
}


#endif RungeKutta

// Do not change the comments below - they will be added automatically by CVS
/*****************************************************************************
*       $Log: RungeKutta.h,v $
*       Revision 1.3  2003/04/23 16:30:59  nilspace
*       Various bugfixes & uploading of all changed code for new programmers.
*       
*       Revision 1.2  2003/03/26 16:38:58  nilspace
*       Updated comments to better work with Doxygen.
*       
*       Revision 1.1  2003/03/05 20:41:04  nilspace
*       Initial submission of example source code file.
*       
*
******************************************************************************/

---