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.math4.analysis.integration.gauss;
018
019import org.apache.commons.math4.analysis.UnivariateFunction;
020import org.apache.commons.math4.exception.DimensionMismatchException;
021import org.apache.commons.math4.util.MathArrays;
022import org.apache.commons.math4.util.Pair;
023
024/**
025 * Class that implements the Gaussian rule for
026 * {@link #integrate(UnivariateFunction) integrating} a weighted
027 * function.
028 *
029 * @since 3.1
030 */
031public class GaussIntegrator {
032    /** Nodes. */
033    private final double[] points;
034    /** Nodes weights. */
035    private final double[] weights;
036
037    /**
038     * Creates an integrator from the given {@code points} and {@code weights}.
039     * The integration interval is defined by the first and last value of
040     * {@code points} which must be sorted in increasing order.
041     *
042     * @param points Integration points.
043     * @param weights Weights of the corresponding integration nodes.
044     * @throws org.apache.commons.math4.exception.NonMonotonicSequenceException if the {@code points} are not
045     * sorted in increasing order.
046     * @throws DimensionMismatchException if points and weights don't have the same length
047     */
048    public GaussIntegrator(double[] points,
049                           double[] weights) {
050        if (points.length != weights.length) {
051            throw new DimensionMismatchException(points.length,
052                                                 weights.length);
053        }
054
055        MathArrays.checkOrder(points, MathArrays.OrderDirection.INCREASING, true, true);
056
057        this.points = points.clone();
058        this.weights = weights.clone();
059    }
060
061    /**
062     * Creates an integrator from the given pair of points (first element of
063     * the pair) and weights (second element of the pair.
064     *
065     * @param pointsAndWeights Integration points and corresponding weights.
066     * @throws org.apache.commons.math4.exception.NonMonotonicSequenceException if the {@code points} are not
067     * sorted in increasing order.
068     *
069     * @see #GaussIntegrator(double[], double[])
070     */
071    public GaussIntegrator(Pair<double[], double[]> pointsAndWeights) {
072        this(pointsAndWeights.getFirst(), pointsAndWeights.getSecond());
073    }
074
075    /**
076     * Returns an estimate of the integral of {@code f(x) * w(x)},
077     * where {@code w} is a weight function that depends on the actual
078     * flavor of the Gauss integration scheme.
079     * The algorithm uses the points and associated weights, as passed
080     * to the {@link #GaussIntegrator(double[],double[]) constructor}.
081     *
082     * @param f Function to integrate.
083     * @return the integral of the weighted function.
084     */
085    public double integrate(UnivariateFunction f) {
086        double s = 0;
087        double c = 0;
088        for (int i = 0; i < points.length; i++) {
089            final double x = points[i];
090            final double w = weights[i];
091            final double y = w * f.value(x) - c;
092            final double t = s + y;
093            c = (t - s) - y;
094            s = t;
095        }
096        return s;
097    }
098
099    /**
100     * @return the order of the integration rule (the number of integration
101     * points).
102     */
103    public int getNumberOfPoints() {
104        return points.length;
105    }
106
107    /**
108     * Gets the integration point at the given index.
109     * The index must be in the valid range but no check is performed.
110     * @param index index of the integration point
111     * @return the integration point.
112     */
113    public double getPoint(int index) {
114        return points[index];
115    }
116
117    /**
118     * Gets the weight of the integration point at the given index.
119     * The index must be in the valid range but no check is performed.
120     * @param index index of the integration point
121     * @return the weight.
122     */
123    public double getWeight(int index) {
124        return weights[index];
125    }
126}