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