GaussIntegrator.java

  1. /*
  2.  * Licensed to the Apache Software Foundation (ASF) under one or more
  3.  * contributor license agreements.  See the NOTICE file distributed with
  4.  * this work for additional information regarding copyright ownership.
  5.  * The ASF licenses this file to You under the Apache License, Version 2.0
  6.  * (the "License"); you may not use this file except in compliance with
  7.  * the License.  You may obtain a copy of the License at
  8.  *
  9.  *      http://www.apache.org/licenses/LICENSE-2.0
  10.  *
  11.  * Unless required by applicable law or agreed to in writing, software
  12.  * distributed under the License is distributed on an "AS IS" BASIS,
  13.  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  14.  * See the License for the specific language governing permissions and
  15.  * limitations under the License.
  16.  */
  17. package org.apache.commons.math4.legacy.analysis.integration.gauss;

  18. import org.apache.commons.math4.legacy.analysis.UnivariateFunction;
  19. import org.apache.commons.math4.legacy.exception.DimensionMismatchException;
  20. import org.apache.commons.math4.legacy.core.MathArrays;
  21. import org.apache.commons.math4.legacy.core.Pair;

  22. /**
  23.  * Class that implements the Gaussian rule for
  24.  * {@link #integrate(UnivariateFunction) integrating} a weighted
  25.  * function.
  26.  *
  27.  * @since 3.1
  28.  */
  29. public class GaussIntegrator {
  30.     /** Nodes. */
  31.     private final double[] points;
  32.     /** Nodes weights. */
  33.     private final double[] weights;

  34.     /**
  35.      * Creates an integrator from the given {@code points} and {@code weights}.
  36.      * The integration interval is defined by the first and last value of
  37.      * {@code points} which must be sorted in increasing order.
  38.      *
  39.      * @param points Integration points.
  40.      * @param weights Weights of the corresponding integration nodes.
  41.      * @throws org.apache.commons.math4.legacy.exception.NonMonotonicSequenceException if the {@code points} are not
  42.      * sorted in increasing order.
  43.      * @throws DimensionMismatchException if points and weights don't have the same length
  44.      */
  45.     public GaussIntegrator(double[] points,
  46.                            double[] weights) {
  47.         if (points.length != weights.length) {
  48.             throw new DimensionMismatchException(points.length,
  49.                                                  weights.length);
  50.         }

  51.         MathArrays.checkOrder(points, MathArrays.OrderDirection.INCREASING, true, true);

  52.         this.points = points.clone();
  53.         this.weights = weights.clone();
  54.     }

  55.     /**
  56.      * Creates an integrator from the given pair of points (first element of
  57.      * the pair) and weights (second element of the pair.
  58.      *
  59.      * @param pointsAndWeights Integration points and corresponding weights.
  60.      * @throws org.apache.commons.math4.legacy.exception.NonMonotonicSequenceException if the {@code points} are not
  61.      * sorted in increasing order.
  62.      *
  63.      * @see #GaussIntegrator(double[], double[])
  64.      */
  65.     public GaussIntegrator(Pair<double[], double[]> pointsAndWeights) {
  66.         this(pointsAndWeights.getFirst(), pointsAndWeights.getSecond());
  67.     }

  68.     /**
  69.      * Returns an estimate of the integral of {@code f(x) * w(x)},
  70.      * where {@code w} is a weight function that depends on the actual
  71.      * flavor of the Gauss integration scheme.
  72.      * The algorithm uses the points and associated weights, as passed
  73.      * to the {@link #GaussIntegrator(double[],double[]) constructor}.
  74.      *
  75.      * @param f Function to integrate.
  76.      * @return the integral of the weighted function.
  77.      */
  78.     public double integrate(UnivariateFunction f) {
  79.         double s = 0;
  80.         double c = 0;
  81.         for (int i = 0; i < points.length; i++) {
  82.             final double x = points[i];
  83.             final double w = weights[i];
  84.             final double y = w * f.value(x) - c;
  85.             final double t = s + y;
  86.             c = (t - s) - y;
  87.             s = t;
  88.         }
  89.         return s;
  90.     }

  91.     /**
  92.      * @return the order of the integration rule (the number of integration
  93.      * points).
  94.      */
  95.     public int getNumberOfPoints() {
  96.         return points.length;
  97.     }

  98.     /**
  99.      * Gets the integration point at the given index.
  100.      * The index must be in the valid range but no check is performed.
  101.      * @param index index of the integration point
  102.      * @return the integration point.
  103.      */
  104.     public double getPoint(int index) {
  105.         return points[index];
  106.     }

  107.     /**
  108.      * Gets the weight of 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 weight.
  112.      */
  113.     public double getWeight(int index) {
  114.         return weights[index];
  115.     }
  116. }