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.dfp.DfpField;
25  import org.hipparchus.exception.MathIllegalArgumentException;
26  import org.hipparchus.util.Pair;
27  
28  /**
29   * Class that provides different ways to compute the nodes and weights to be
30   * used by the {@link GaussIntegrator Gaussian integration rule}.
31   */
32  public class GaussIntegratorFactory {
33  
34      /** Number of digits for Legendre high precision. */
35      public static final int DEFAULT_DECIMAL_DIGITS = 40;
36  
37      /** Generator of Gauss-Legendre integrators. */
38      private final RuleFactory legendre;
39      /** Generator of Gauss-Legendre integrators. */
40      private final RuleFactory legendreHighPrecision;
41      /** Generator of Gauss-Hermite integrators. */
42      private final RuleFactory hermite;
43      /** Generator of Gauss-Laguerre integrators. */
44      private final RuleFactory laguerre;
45  
46      /** Simple constructor.
47       */
48      public GaussIntegratorFactory() {
49          this(DEFAULT_DECIMAL_DIGITS);
50      }
51  
52      /** Simple constructor.
53       * @param decimalDigits minimum number of decimal digits for {@link #legendreHighPrecision(int)}
54       */
55      public GaussIntegratorFactory(final int decimalDigits) {
56          legendre              = new LegendreRuleFactory();
57          legendreHighPrecision = new ConvertingRuleFactory<>(new FieldLegendreRuleFactory<>(new DfpField(decimalDigits)));
58          hermite               = new HermiteRuleFactory();
59          laguerre              = new LaguerreRuleFactory();
60      }
61  
62      /**
63       * Creates a Gauss-Laguerre integrator of the given order.
64       * The call to the
65       * {@link GaussIntegrator#integrate(org.hipparchus.analysis.UnivariateFunction)
66       * integrate} method will perform an integration on the interval
67       * \([0, +\infty)\): the computed value is the improper integral of
68       * \(e^{-x} f(x)\)
69       * where \(f(x)\) is the function passed to the
70       * {@link SymmetricGaussIntegrator#integrate(org.hipparchus.analysis.UnivariateFunction)
71       * integrate} method.
72       *
73       * @param numberOfPoints Order of the integration rule.
74       * @return a Gauss-Legendre integrator.
75       */
76      public GaussIntegrator laguerre(int numberOfPoints) {
77          return new GaussIntegrator(laguerre.getRule(numberOfPoints));
78      }
79  
80      /**
81       * Creates a Gauss-Legendre integrator of the given order.
82       * The call to the
83       * {@link GaussIntegrator#integrate(org.hipparchus.analysis.UnivariateFunction)
84       * integrate} method will perform an integration on the natural interval
85       * {@code [-1 , 1]}.
86       *
87       * @param numberOfPoints Order of the integration rule.
88       * @return a Gauss-Legendre integrator.
89       */
90      public GaussIntegrator legendre(int numberOfPoints) {
91          return new GaussIntegrator(legendre.getRule(numberOfPoints));
92      }
93  
94      /**
95       * Creates a Gauss-Legendre integrator of the given order.
96       * The call to the
97       * {@link GaussIntegrator#integrate(org.hipparchus.analysis.UnivariateFunction)
98       * integrate} method will perform an integration on the given interval.
99       *
100      * @param numberOfPoints Order of the integration rule.
101      * @param lowerBound Lower bound of the integration interval.
102      * @param upperBound Upper bound of the integration interval.
103      * @return a Gauss-Legendre integrator.
104      * @throws MathIllegalArgumentException if number of points is not positive
105      */
106     public GaussIntegrator legendre(int numberOfPoints,
107                                     double lowerBound,
108                                     double upperBound)
109         throws MathIllegalArgumentException {
110         return new GaussIntegrator(transform(legendre.getRule(numberOfPoints),
111                                              lowerBound, upperBound));
112     }
113 
114     /**
115      * Creates a Gauss-Legendre integrator of the given order.
116      * The call to the
117      * {@link GaussIntegrator#integrate(org.hipparchus.analysis.UnivariateFunction)
118      * integrate} method will perform an integration on the natural interval
119      * {@code [-1 , 1]}.
120      *
121      * @param numberOfPoints Order of the integration rule.
122      * @return a Gauss-Legendre integrator.
123      * @throws MathIllegalArgumentException if number of points is not positive
124      */
125     public GaussIntegrator legendreHighPrecision(int numberOfPoints)
126         throws MathIllegalArgumentException {
127         return new GaussIntegrator(legendreHighPrecision.getRule(numberOfPoints));
128     }
129 
130     /**
131      * Creates an integrator of the given order, and whose call to the
132      * {@link GaussIntegrator#integrate(org.hipparchus.analysis.UnivariateFunction)
133      * integrate} method will perform an integration on the given interval.
134      *
135      * @param numberOfPoints Order of the integration rule.
136      * @param lowerBound Lower bound of the integration interval.
137      * @param upperBound Upper bound of the integration interval.
138      * @return a Gauss-Legendre integrator.
139      * @throws MathIllegalArgumentException if number of points is not positive
140      */
141     public GaussIntegrator legendreHighPrecision(int numberOfPoints,
142                                                  double lowerBound,
143                                                  double upperBound)
144         throws MathIllegalArgumentException {
145         return new GaussIntegrator(transform(legendreHighPrecision.getRule(numberOfPoints),
146                                              lowerBound, upperBound));
147     }
148 
149     /**
150      * Creates a Gauss-Hermite integrator of the given order.
151      * The call to the
152      * {@link SymmetricGaussIntegrator#integrate(org.hipparchus.analysis.UnivariateFunction)
153      * integrate} method will perform a weighted integration on the interval
154      * \([-\infty, +\infty]\): the computed value is the improper integral of
155      * \(e^{-x^2}f(x)\)
156      * where \(f(x)\) is the function passed to the
157      * {@link SymmetricGaussIntegrator#integrate(org.hipparchus.analysis.UnivariateFunction)
158      * integrate} method.
159      *
160      * @param numberOfPoints Order of the integration rule.
161      * @return a Gauss-Hermite integrator.
162      */
163     public SymmetricGaussIntegrator hermite(int numberOfPoints) {
164         return new SymmetricGaussIntegrator(hermite.getRule(numberOfPoints));
165     }
166 
167     /**
168      * Performs a change of variable so that the integration can be performed
169      * on an arbitrary interval {@code [a, b]}.
170      * It is assumed that the natural interval is {@code [-1, 1]}.
171      *
172      * @param rule Original points and weights.
173      * @param a Lower bound of the integration interval.
174      * @param b Lower bound of the integration interval.
175      * @return the points and weights adapted to the new interval.
176      */
177     private Pair<double[], double[]> transform(Pair<double[], double[]> rule, double a, double b) {
178         final double[] points = rule.getFirst();
179         final double[] weights = rule.getSecond();
180 
181         // Scaling
182         final double scale = (b - a) / 2;
183         final double shift = a + scale;
184 
185         for (int i = 0; i < points.length; i++) {
186             points[i] = points[i] * scale + shift;
187             weights[i] *= scale;
188         }
189 
190         return new Pair<>(points, weights);
191     }
192 
193 }