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.CalculusFieldElement;
25 import org.hipparchus.analysis.CalculusFieldUnivariateFunction;
26 import org.hipparchus.exception.MathIllegalArgumentException;
27 import org.hipparchus.util.Pair;
28
29 /**
30 * This class's implements {@link #integrate(CalculusFieldUnivariateFunction) integrate}
31 * method assuming that the integral is symmetric about 0.
32 * This allows to reduce numerical errors.
33 *
34 * @param <T> Type of the field elements.
35 * @since 2.0
36 */
37 public class SymmetricFieldGaussIntegrator<T extends CalculusFieldElement<T>> extends FieldGaussIntegrator<T> {
38 /**
39 * Creates an integrator from the given {@code points} and {@code weights}.
40 * The integration interval is defined by the first and last value of
41 * {@code points} which must be sorted in increasing order.
42 *
43 * @param points Integration points.
44 * @param weights Weights of the corresponding integration nodes.
45 * @throws MathIllegalArgumentException if the {@code points} are not
46 * sorted in increasing order.
47 * @throws MathIllegalArgumentException if points and weights don't have the same length
48 */
49 public SymmetricFieldGaussIntegrator(T[] points, T[] weights)
50 throws MathIllegalArgumentException {
51 super(points, weights);
52 }
53
54 /**
55 * Creates an integrator from the given pair of points (first element of
56 * the pair) and weights (second element of the pair.
57 *
58 * @param pointsAndWeights Integration points and corresponding weights.
59 * @throws MathIllegalArgumentException if the {@code points} are not
60 * sorted in increasing order.
61 *
62 * @see #SymmetricFieldGaussIntegrator(CalculusFieldElement[], CalculusFieldElement[])
63 */
64 public SymmetricFieldGaussIntegrator(Pair<T[], T[]> pointsAndWeights)
65 throws MathIllegalArgumentException {
66 this(pointsAndWeights.getFirst(), pointsAndWeights.getSecond());
67 }
68
69 /**
70 * {@inheritDoc}
71 */
72 @Override
73 public T integrate(CalculusFieldUnivariateFunction<T> f) {
74 final int ruleLength = getNumberOfPoints();
75
76 final T zero = getPoint(0).getField().getZero();
77 if (ruleLength == 1) {
78 return getWeight(0).multiply(f.value(zero));
79 }
80
81 final int iMax = ruleLength / 2;
82 T s = zero;
83 T c = zero;
84 for (int i = 0; i < iMax; i++) {
85 final T p = getPoint(i);
86 final T w = getWeight(i);
87
88 final T f1 = f.value(p);
89 final T f2 = f.value(p.negate());
90
91 final T y = w.multiply(f1.add(f2)).subtract(c);
92 final T t = s.add(y);
93
94 c = t.subtract(s).subtract(y);
95 s = t;
96 }
97
98 if (ruleLength % 2 != 0) {
99 final T w = getWeight(iMax);
100
101 final T y = w.multiply(f.value(zero)).subtract(c);
102 final T t = s.add(y);
103
104 s = t;
105 }
106
107 return s;
108 }
109 }