RungeKuttaIntegrator1D.java

/**
 * Copyright (C) 2009 - present by OpenGamma Inc. and the OpenGamma group of companies
 *
 * Please see distribution for license.
 */
package com.opengamma.analytics.math.integration;

import com.opengamma.OpenGammaRuntimeException;
import com.opengamma.analytics.math.function.Function1D;
import com.opengamma.util.ArgumentChecker;

/**
 * Adapted from the forth-order Runge-Kutta method for solving ODE. See <a
 * href="http://en.wikipedia.org/wiki/Runge-Kutta_methods">here </a> for the
 * maths. It is a very robust integrator and should be used before trying more
 * specialised methods.
 */
public class RungeKuttaIntegrator1D extends Integrator1D<Double, Double> {
  private static final double DEF_TOL = 1e-10;
  private static final double STEP_SIZE_LIMIT = 1e-50;
  private static final int DEF_MIN_STEPS = 10;
  private final double _absTol, _relTol;
  private final int _minSteps;

  public RungeKuttaIntegrator1D(final double absTol, final double relTol, final int minSteps) {
    if (absTol < 0.0 || Double.isNaN(absTol) || Double.isInfinite(absTol)) {
      throw new IllegalArgumentException("Absolute Tolerance must be greater than zero");
    }
    if (relTol < 0.0 || Double.isNaN(relTol) || Double.isInfinite(relTol)) {
      throw new IllegalArgumentException("Relative Tolerance must be greater than zero");
    }
    if (minSteps < 1) {
      throw new IllegalArgumentException("Must have minimum of 1 step");
    }
    _absTol = absTol;
    _relTol = relTol;
    _minSteps = minSteps;
  }

  public RungeKuttaIntegrator1D(final double tol, final int minSteps) {
    this(tol, tol, minSteps);
  }

  public RungeKuttaIntegrator1D(final double atol, final double rtol) {
    this(atol, rtol, DEF_MIN_STEPS);
  }

  public RungeKuttaIntegrator1D(final double tol) {
    this(tol, tol, DEF_MIN_STEPS);
  }

  public RungeKuttaIntegrator1D(final int minSteps) {
    this(DEF_TOL, minSteps);
  }

  public RungeKuttaIntegrator1D() {
    this(DEF_TOL, DEF_MIN_STEPS);

  }

  @Override
  public Double integrate(final Function1D<Double, Double> f, final Double lower, final Double upper) {
    ArgumentChecker.notNull(lower, "lower");
    ArgumentChecker.notNull(upper, "upper");
    if (Double.isNaN(lower) || Double.isInfinite(lower) || Double.isInfinite(upper) || Double.isNaN(upper)) {
      throw new OpenGammaRuntimeException("lower or upper was NaN or Inf");
    }
    
    final double h = (upper - lower) / _minSteps;
    double f1, f2, f3, x;
    x = lower;
    f1 = f.evaluate(x);
    if (Double.isNaN(f1) || Double.isInfinite(f1)) {
      throw new OpenGammaRuntimeException("function evaluation returned NaN or Inf");
    }

    double result = 0.0;
    for (int i = 0; i < _minSteps; i++) {
      f2 = f.evaluate(x + h / 2.0);
      if (Double.isNaN(f2) || Double.isInfinite(f2)) {
        throw new OpenGammaRuntimeException("function evaluation returned NaN or Inf");
      }
      f3 = f.evaluate(x + h);
      if (Double.isNaN(f3) || Double.isInfinite(f3)) {
        throw new OpenGammaRuntimeException("function evaluation returned NaN or Inf");
      }

      result += calculateRungeKuttaFourthOrder(f, x, h, f1, f2, f3);
      f1 = f3;
      x += h;
    }
    return result;
  }

  private double calculateRungeKuttaFourthOrder(final Function1D<Double, Double> f, final double x, final double h, final double fl, final double fm, final double fu) {
//    if (Double.isNaN(h) || Double.isInfinite(h) || 
//        Double.isNaN(fl) || Double.isInfinite(fl) ||
//        Double.isNaN(fm) || Double.isInfinite(fm) ||
//        Double.isNaN(fu) || Double.isInfinite(fu)) {
//      throw new OpenGammaRuntimeException("h was Inf or NaN");
//    }
    final double f1 = f.evaluate(x + 0.25 * h);
    if (Double.isNaN(f1) || Double.isInfinite(f1)) {
      throw new OpenGammaRuntimeException("f.evaluate returned NaN or Inf");
    }
    final double f2 = f.evaluate(x + 0.75 * h);
    if (Double.isNaN(f2) || Double.isInfinite(f2)) {
      throw new OpenGammaRuntimeException("f.evaluate returned NaN or Inf");
    }
    final double ya = h * (fl + 4.0 * fm + fu) / 6.0;
    final double yb = h * (fl + 2.0 * fm + 4.0 * (f1 + f2) + fu) / 12.0;

    final double diff = Math.abs(ya - yb);
    final double abs = Math.max(Math.abs(ya), Math.abs(yb));

    if (diff < _absTol + _relTol * abs) {
      return yb + (yb - ya) / 15.0;
    }

    // can't keep halving the step size
    if (h < STEP_SIZE_LIMIT) {
      return yb + (yb - ya) / 15.0;
    }

    return calculateRungeKuttaFourthOrder(f, x, h / 2.0, fl, f1, fm) + calculateRungeKuttaFourthOrder(f, x + h / 2.0, h / 2.0, fm, f2, fu);
  }

}