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