/**
* 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 org.apache.commons.lang.Validate;
import org.apache.commons.math.util.MathUtils;
import com.opengamma.analytics.math.function.DoubleFunction1D;
import com.opengamma.analytics.math.function.Function1D;
import com.opengamma.analytics.math.function.special.GammaFunction;
import com.opengamma.analytics.math.function.special.LaguerrePolynomialFunction;
import com.opengamma.analytics.math.rootfinding.NewtonRaphsonSingleRootFinder;
import com.opengamma.util.tuple.Pair;
/**
* Class that generates weights and abscissas for Gauss-Laguerre quadrature.
* The weights $w_i$ are given by:
* $$
* \begin{align*}
* w_i = -\frac{\Gamma(\alpha + n)}{n!L_i'(x_i)L_{i-1}(x_i)}
* \end{align*}
* $$
* where $x_i$ is the $i^{th}$ root of the orthogonal polynomial, $L_i$ is the
* $i^{th}$ polynomial and $L_i'$ is the first derivative of the $i^{th}$
* polynomial. The orthogonal polynomial is generated by
* {@link com.opengamma.analytics.math.function.special.LaguerrePolynomialFunction}.
*/
public class GaussLaguerreWeightAndAbscissaFunction implements QuadratureWeightAndAbscissaFunction {
private static final LaguerrePolynomialFunction LAGUERRE = new LaguerrePolynomialFunction();
private static final NewtonRaphsonSingleRootFinder ROOT_FINDER = new NewtonRaphsonSingleRootFinder(1e-10);
private static final Function1D<Double, Double> GAMMA_FUNCTION = new GammaFunction();
private final double _alpha;
/**
* Sets $\alpha = 0$
*/
public GaussLaguerreWeightAndAbscissaFunction() {
this(0);
}
/**
* @param alpha The value of $\alpha$ to use when generating the polynomials.
*/
public GaussLaguerreWeightAndAbscissaFunction(final double alpha) {
_alpha = alpha;
}
/**
* {@inheritDoc}
*/
@Override
public GaussianQuadratureData generate(final int n) {
Validate.isTrue(n > 0);
final Pair<DoubleFunction1D, DoubleFunction1D>[] polynomials = LAGUERRE.getPolynomialsAndFirstDerivative(n, _alpha);
final Pair<DoubleFunction1D, DoubleFunction1D> pair = polynomials[n];
final DoubleFunction1D p1 = polynomials[n - 1].getFirst();
final DoubleFunction1D function = pair.getFirst();
final DoubleFunction1D derivative = pair.getSecond();
final double[] x = new double[n];
final double[] w = new double[n];
double root = 0;
for (int i = 0; i < n; i++) {
root = ROOT_FINDER.getRoot(function, derivative, getInitialRootGuess(root, i, n, x));
x[i] = root;
w[i] = -GAMMA_FUNCTION.evaluate(_alpha + n) / MathUtils.factorialDouble(n) / (derivative.evaluate(root) * p1.evaluate(root));
}
return new GaussianQuadratureData(x, w);
}
private double getInitialRootGuess(final double previousRoot, final int i, final int n, final double[] x) {
if (i == 0) {
return (1 + _alpha) * (3 + 0.92 * _alpha) / (1 + 1.8 * _alpha + 2.4 * n);
}
if (i == 1) {
return previousRoot + (15 + 6.25 * _alpha) / (1 + 0.9 * _alpha + 2.5 * n);
}
final int j = i - 1;
return previousRoot + ((1 + 2.55 * j) / 1.9 / j + 1.26 * j * _alpha / (1 + 3.5 * j)) * (previousRoot - x[i - 2]) / (1 + 0.3 * _alpha);
}
}