View Javadoc
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    *      https://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  
18  /*
19   * This is not the original file distributed by the Apache Software Foundation
20   * It has been modified by the Hipparchus project
21   */
22  package org.hipparchus.analysis.integration.gauss;
23  
24  import org.hipparchus.analysis.UnivariateFunction;
25  import org.hipparchus.exception.MathIllegalArgumentException;
26  import org.hipparchus.util.MathArrays;
27  import org.hipparchus.util.Pair;
28  
29  /**
30   * Class that implements the Gaussian rule for
31   * {@link #integrate(UnivariateFunction) integrating} a weighted
32   * function.
33   *
34   */
35  public class GaussIntegrator {
36      /** Nodes. */
37      private final double[] points;
38      /** Nodes weights. */
39      private final double[] weights;
40  
41      /**
42       * Creates an integrator from the given {@code points} and {@code weights}.
43       * The integration interval is defined by the first and last value of
44       * {@code points} which must be sorted in increasing order.
45       *
46       * @param points Integration points.
47       * @param weights Weights of the corresponding integration nodes.
48       * @throws MathIllegalArgumentException if the {@code points} are not
49       * sorted in increasing order.
50       * @throws MathIllegalArgumentException if points and weights don't have the same length
51       */
52      public GaussIntegrator(double[] points,
53                             double[] weights)
54          throws MathIllegalArgumentException {
55  
56          MathArrays.checkEqualLength(points, weights);
57          MathArrays.checkOrder(points, MathArrays.OrderDirection.INCREASING, true, true);
58  
59          this.points = points.clone();
60          this.weights = weights.clone();
61      }
62  
63      /**
64       * Creates an integrator from the given pair of points (first element of
65       * the pair) and weights (second element of the pair.
66       *
67       * @param pointsAndWeights Integration points and corresponding weights.
68       * @throws MathIllegalArgumentException if the {@code points} are not
69       * sorted in increasing order.
70       *
71       * @see #GaussIntegrator(double[], double[])
72       */
73      public GaussIntegrator(Pair<double[], double[]> pointsAndWeights)
74          throws MathIllegalArgumentException {
75          this(pointsAndWeights.getFirst(), pointsAndWeights.getSecond());
76      }
77  
78      /**
79       * Returns an estimate of the integral of {@code f(x) * w(x)},
80       * where {@code w} is a weight function that depends on the actual
81       * flavor of the Gauss integration scheme.
82       * The algorithm uses the points and associated weights, as passed
83       * to the {@link #GaussIntegrator(double[],double[]) constructor}.
84       *
85       * @param f Function to integrate.
86       * @return the integral of the weighted function.
87       */
88      public double integrate(UnivariateFunction f) {
89          double s = 0;
90          double c = 0;
91          for (int i = 0; i < points.length; i++) {
92              final double x = points[i];
93              final double w = weights[i];
94              final double y = w * f.value(x) - c;
95              final double t = s + y;
96              c = (t - s) - y;
97              s = t;
98          }
99          return s;
100     }
101 
102     /** Get the order of the integration rule.
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 }