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.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 }