adaptiverungekutta.cpp

dysii是一款非常出色的滤波函数库

//#if defined(__GNUC__) && defined(GCC_PCH)
//  #include "../aux/aux.hpp"
//#endif

#include "AdaptiveRungeKutta.hpp"

#include <assert.h>
#include <math.h>

#include <gsl/gsl_errno.h>

using namespace indii::ml::ode;

namespace aux = indii::ml::aux;

const gsl_odeiv_step_type* AdaptiveRungeKutta::gslStepType
    = gsl_odeiv_step_rkf45;  // explicit method
  //= gsl_odeiv_step_rk4imp; // implicit method

AdaptiveRungeKutta::AdaptiveRungeKutta(DifferentialModel* model) :
  NumericalSolver(model->getDimensions()), model(model) {
  setForwardFunction(gslForwardFunction);
  setBackwardFunction(gslBackwardFunction);
  setStepType(gslStepType);
}

AdaptiveRungeKutta::AdaptiveRungeKutta(DifferentialModel* model,
    const aux::vector& y0) : NumericalSolver(y0), model(model) {
  /* pre-condition */
  assert (y0.size() == model->getDimensions());

  setForwardFunction(gslForwardFunction);
  setBackwardFunction(gslBackwardFunction);
  setStepType(gslStepType);
}

AdaptiveRungeKutta::~AdaptiveRungeKutta() {
  //
}

int AdaptiveRungeKutta::gslForwardFunction(double t, const double y[],
    double dydt[], void *params) {
  AdaptiveRungeKutta* m = (AdaptiveRungeKutta*)params;
  m->model->calculateDerivatives(t, y, dydt);

  return GSL_SUCCESS;
}

int AdaptiveRungeKutta::gslBackwardFunction(double t, const double y[],
    double dydt[], void *params) {
  AdaptiveRungeKutta* m = (AdaptiveRungeKutta*)params;

  /* convert the proposed step to a future time into a proposed step
     to a past time */
  int result = gslForwardFunction(2.0 * m->base - t, y, dydt, params);

  if (result == GSL_SUCCESS) {
    /* as we're moving backward in time, negate gradients */
    size_t i;
    for (i = 0; i < m->dimensions; i++) {
      dydt[i] *= -1.0;
    }
  }

  return result;
}