    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    }