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