001/*
002 * Licensed to the Apache Software Foundation (ASF) under one or more
003 * contributor license agreements.  See the NOTICE file distributed with
004 * this work for additional information regarding copyright ownership.
005 * The ASF licenses this file to You under the Apache License, Version 2.0
006 * (the "License"); you may not use this file except in compliance with
007 * the License.  You may obtain a copy of the License at
008 *
009 *      http://www.apache.org/licenses/LICENSE-2.0
010 *
011 * Unless required by applicable law or agreed to in writing, software
012 * distributed under the License is distributed on an "AS IS" BASIS,
013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014 * See the License for the specific language governing permissions and
015 * limitations under the License.
016 */
017package org.apache.commons.math3.analysis.integration.gauss;
018
019import org.apache.commons.math3.analysis.UnivariateFunction;
020import org.apache.commons.math3.exception.DimensionMismatchException;
021import org.apache.commons.math3.exception.NonMonotonicSequenceException;
022import org.apache.commons.math3.util.MathArrays;
023import org.apache.commons.math3.util.Pair;
024
025/**
026 * Class that implements the Gaussian rule for
027 * {@link #integrate(UnivariateFunction) integrating} a weighted
028 * function.
029 *
030 * @since 3.1
031 */
032public class GaussIntegrator {
033    /** Nodes. */
034    private final double[] points;
035    /** Nodes weights. */
036    private final double[] weights;
037
038    /**
039     * Creates an integrator from the given {@code points} and {@code weights}.
040     * The integration interval is defined by the first and last value of
041     * {@code points} which must be sorted in increasing order.
042     *
043     * @param points Integration points.
044     * @param weights Weights of the corresponding integration nodes.
045     * @throws NonMonotonicSequenceException if the {@code points} are not
046     * sorted in increasing order.
047     * @throws DimensionMismatchException if points and weights don't have the same length
048     */
049    public GaussIntegrator(double[] points,
050                           double[] weights)
051        throws NonMonotonicSequenceException, DimensionMismatchException {
052        if (points.length != weights.length) {
053            throw new DimensionMismatchException(points.length,
054                                                 weights.length);
055        }
056
057        MathArrays.checkOrder(points, MathArrays.OrderDirection.INCREASING, true, true);
058
059        this.points = points.clone();
060        this.weights = weights.clone();
061    }
062
063    /**
064     * Creates an integrator from the given pair of points (first element of
065     * the pair) and weights (second element of the pair.
066     *
067     * @param pointsAndWeights Integration points and corresponding weights.
068     * @throws NonMonotonicSequenceException if the {@code points} are not
069     * sorted in increasing order.
070     *
071     * @see #GaussIntegrator(double[], double[])
072     */
073    public GaussIntegrator(Pair<double[], double[]> pointsAndWeights)
074        throws NonMonotonicSequenceException {
075        this(pointsAndWeights.getFirst(), pointsAndWeights.getSecond());
076    }
077
078    /**
079     * Returns an estimate of the integral of {@code f(x) * w(x)},
080     * where {@code w} is a weight function that depends on the actual
081     * flavor of the Gauss integration scheme.
082     * The algorithm uses the points and associated weights, as passed
083     * to the {@link #GaussIntegrator(double[],double[]) constructor}.
084     *
085     * @param f Function to integrate.
086     * @return the integral of the weighted function.
087     */
088    public double integrate(UnivariateFunction f) {
089        double s = 0;
090        double c = 0;
091        for (int i = 0; i < points.length; i++) {
092            final double x = points[i];
093            final double w = weights[i];
094            final double y = w * f.value(x) - c;
095            final double t = s + y;
096            c = (t - s) - y;
097            s = t;
098        }
099        return s;
100    }
101
102    /**
103     * @return the order of the integration rule (the number of integration
104     * points).
105     */
106    public int getNumberOfPoints() {
107        return points.length;
108    }
109
110    /**
111     * Gets the integration point at the given index.
112     * The index must be in the valid range but no check is performed.
113     * @param index index of the integration point
114     * @return the integration point.
115     */
116    public double getPoint(int index) {
117        return points[index];
118    }
119
120    /**
121     * Gets the weight of the integration point at the given index.
122     * The index must be in the valid range but no check is performed.
123     * @param index index of the integration point
124     * @return the weight.
125     */
126    public double getWeight(int index) {
127        return weights[index];
128    }
129}