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 }