RungeKutta4Solver.cpp

/***************************************************************************
 *   This file is part of OpenCAL: Open Computer Animation Library         *
 *   I created OpenCAL as my master's thesis Computer Science (multimedia) *
 *   at the tUL university in Diepenbeek, Belgium                          *
 *                                                                         *
 *   Copyright (C) 2003-2004 by Jeroen Dierckx                             *
 *   jeroen.dierckx@student.luc.ac.be                                      *
 *                                                                         *
 ***************************************************************************/

// Includes
#include "RungeKutta4Solver.h"
#include <OpenCAL/ODESource.h>
using namespace OpenCAL;

using namespace std;


/******************************
* Constructors and destructor *
******************************/

RungeKutta4Solver::RungeKutta4Solver()
{
}

RungeKutta4Solver::~RungeKutta4Solver()
{
}


/************************
* Get and set functions *
************************/


/******************
* Other functions *
******************/

void RungeKutta4Solver::solve(float stepSize, ODESource *source)
{
    const float oneThird = 1.0f / 3.0f;
    const float oneSixth = 1.0f / 6.0f;

    source->fillState();
    source->fillDerivative();

    PhaseSpace &state = *(source->getState());
    PhaseSpace &deriv = *(source->getDerivative());

    // Copy the original state
    PhaseSpace x0 = state;


    // Calculate k1
    PhaseSpace k1 = stepSize * deriv;

    // Calculate k2
    state.add(k1, 0.5f);
    source->fromState();
    source->fillDerivative();
    PhaseSpace k2 = stepSize * deriv;

    // Calculate k3
    state = x0;
    state.add(k2, 0.5f);
    source->fromState();
    source->fillDerivative();
    PhaseSpace k3 = stepSize * deriv;

    // Calculate k4
    state = x0;
    state += k3;
    source->fromState();
    source->fillDerivative();
    PhaseSpace k4 = stepSize * deriv;


    // Now calculate the new state

    // This works, but requires 9 copies of a phase space
    //state = x0 + oneSixth * k1 + oneThird * k2 + oneThird * k3 + oneSixth * k4;

    // This requires only one copy
    state = x0;
    state.add(k1, oneSixth);
    state.add(k2, oneThird);
    state.add(k3, oneThird);
    state.add(k4, oneSixth);


    source->fromState();
}

---