com.opengamma.analytics.math.integration.GaussLegendreWeightAndAbscissaFunction.java Source code

Introduction

Here is the source code for com.opengamma.analytics.math.integration.GaussLegendreWeightAndAbscissaFunction.java

Source

/**
 * 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));
    }
}