/**
* 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 com.opengamma.analytics.math.function.DoubleFunction1D;
import com.opengamma.analytics.math.function.special.LegendrePolynomialFunction;
import com.opengamma.analytics.math.rootfinding.NewtonRaphsonSingleRootFinder;
import com.opengamma.util.tuple.Pair;
/**
* Class that generates weights and abscissas for Gauss-Legendre quadrature.
* The weights $w_i$ are given by:
* $$
* \begin{align*}
* w_i = \frac{2}{(1 - x_i^2) L_i'(x_i)^2}
* \end{align*}
* $$
* where $x_i$ is the $i^{th}$ root of the orthogonal 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.LegendrePolynomialFunction}.
*/
public class GaussLegendreWeightAndAbscissaFunction implements QuadratureWeightAndAbscissaFunction {
private static final LegendrePolynomialFunction LEGENDRE = new LegendrePolynomialFunction();
private static final NewtonRaphsonSingleRootFinder ROOT_FINDER = new NewtonRaphsonSingleRootFinder(1e-15);
/**
* {@inheritDoc}
*/
@Override
public GaussianQuadratureData generate(final int n) {
Validate.isTrue(n > 0);
final int mid = (n + 1) / 2;
final double[] x = new double[n];
final double[] w = new double[n];
final Pair<DoubleFunction1D, DoubleFunction1D>[] polynomials = LEGENDRE.getPolynomialsAndFirstDerivative(n);
final Pair<DoubleFunction1D, DoubleFunction1D> pair = polynomials[n];
final DoubleFunction1D function = pair.getFirst();
final DoubleFunction1D derivative = pair.getSecond();
for (int i = 0; i < mid; i++) {
final double root = ROOT_FINDER.getRoot(function, derivative, getInitialRootGuess(i, n));
x[i] = -root;
x[n - i - 1] = root;
final double dp = derivative.evaluate(root);
w[i] = 2 / ((1 - root * root) * dp * dp);
w[n - i - 1] = w[i];
}
return new GaussianQuadratureData(x, w);
}
private double getInitialRootGuess(final int i, final int n) {
return Math.cos(Math.PI * (i + 0.75) / (n + 0.5));
}
}