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     */
017    package org.apache.commons.math3.analysis.integration.gauss;
018    
019    import org.apache.commons.math3.analysis.UnivariateFunction;
020    import org.apache.commons.math3.exception.DimensionMismatchException;
021    import org.apache.commons.math3.exception.NonMonotonicSequenceException;
022    import org.apache.commons.math3.util.MathArrays;
023    import 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     * @version $Id: GaussIntegrator.java 1382197 2012-09-07 22:35:01Z erans $
032     */
033    public class GaussIntegrator {
034        /** Nodes. */
035        private final double[] points;
036        /** Nodes weights. */
037        private final double[] weights;
038    
039        /**
040         * Creates an integrator from the given {@code points} and {@code weights}.
041         * The integration interval is defined by the first and last value of
042         * {@code points} which must be sorted in increasing order.
043         *
044         * @param points Integration points.
045         * @param weights Weights of the corresponding integration nodes.
046         * @throws NonMonotonicSequenceException if the {@code points} are not
047         * sorted in increasing order.
048         */
049        public GaussIntegrator(double[] points,
050                               double[] weights)
051            throws NonMonotonicSequenceException {
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    }