1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
|
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math3.analysis.integration.gauss;
import java.math.BigDecimal;
import org.apache.commons.math3.exception.DimensionMismatchException;
import org.apache.commons.math3.exception.NotStrictlyPositiveException;
import org.apache.commons.math3.util.Pair;
/**
* Class that provides different ways to compute the nodes and weights to be
* used by the {@link GaussIntegrator Gaussian integration rule}.
*
* @since 3.1
*/
public class GaussIntegratorFactory {
/** Generator of Gauss-Legendre integrators. */
private final BaseRuleFactory<Double> legendre = new LegendreRuleFactory();
/** Generator of Gauss-Legendre integrators. */
private final BaseRuleFactory<BigDecimal> legendreHighPrecision = new LegendreHighPrecisionRuleFactory();
/** Generator of Gauss-Hermite integrators. */
private final BaseRuleFactory<Double> hermite = new HermiteRuleFactory();
/**
* Creates a Gauss-Legendre integrator of the given order.
* The call to the
* {@link GaussIntegrator#integrate(org.apache.commons.math3.analysis.UnivariateFunction)
* integrate} method will perform an integration on the natural interval
* {@code [-1 , 1]}.
*
* @param numberOfPoints Order of the integration rule.
* @return a Gauss-Legendre integrator.
*/
public GaussIntegrator legendre(int numberOfPoints) {
return new GaussIntegrator(getRule(legendre, numberOfPoints));
}
/**
* Creates a Gauss-Legendre integrator of the given order.
* The call to the
* {@link GaussIntegrator#integrate(org.apache.commons.math3.analysis.UnivariateFunction)
* integrate} method will perform an integration on the given interval.
*
* @param numberOfPoints Order of the integration rule.
* @param lowerBound Lower bound of the integration interval.
* @param upperBound Upper bound of the integration interval.
* @return a Gauss-Legendre integrator.
* @throws NotStrictlyPositiveException if number of points is not positive
*/
public GaussIntegrator legendre(int numberOfPoints,
double lowerBound,
double upperBound)
throws NotStrictlyPositiveException {
return new GaussIntegrator(transform(getRule(legendre, numberOfPoints),
lowerBound, upperBound));
}
/**
* Creates a Gauss-Legendre integrator of the given order.
* The call to the
* {@link GaussIntegrator#integrate(org.apache.commons.math3.analysis.UnivariateFunction)
* integrate} method will perform an integration on the natural interval
* {@code [-1 , 1]}.
*
* @param numberOfPoints Order of the integration rule.
* @return a Gauss-Legendre integrator.
* @throws NotStrictlyPositiveException if number of points is not positive
*/
public GaussIntegrator legendreHighPrecision(int numberOfPoints)
throws NotStrictlyPositiveException {
return new GaussIntegrator(getRule(legendreHighPrecision, numberOfPoints));
}
/**
* Creates an integrator of the given order, and whose call to the
* {@link GaussIntegrator#integrate(org.apache.commons.math3.analysis.UnivariateFunction)
* integrate} method will perform an integration on the given interval.
*
* @param numberOfPoints Order of the integration rule.
* @param lowerBound Lower bound of the integration interval.
* @param upperBound Upper bound of the integration interval.
* @return a Gauss-Legendre integrator.
* @throws NotStrictlyPositiveException if number of points is not positive
*/
public GaussIntegrator legendreHighPrecision(int numberOfPoints,
double lowerBound,
double upperBound)
throws NotStrictlyPositiveException {
return new GaussIntegrator(transform(getRule(legendreHighPrecision, numberOfPoints),
lowerBound, upperBound));
}
/**
* Creates a Gauss-Hermite integrator of the given order.
* The call to the
* {@link SymmetricGaussIntegrator#integrate(org.apache.commons.math3.analysis.UnivariateFunction)
* integrate} method will perform a weighted integration on the interval
* \([-\infty, +\infty]\): the computed value is the improper integral of
* \(e^{-x^2}f(x)\)
* where \(f(x)\) is the function passed to the
* {@link SymmetricGaussIntegrator#integrate(org.apache.commons.math3.analysis.UnivariateFunction)
* integrate} method.
*
* @param numberOfPoints Order of the integration rule.
* @return a Gauss-Hermite integrator.
*/
public SymmetricGaussIntegrator hermite(int numberOfPoints) {
return new SymmetricGaussIntegrator(getRule(hermite, numberOfPoints));
}
/**
* @param factory Integration rule factory.
* @param numberOfPoints Order of the integration rule.
* @return the integration nodes and weights.
* @throws NotStrictlyPositiveException if number of points is not positive
* @throws DimensionMismatchException if the elements of the rule pair do not
* have the same length.
*/
private static Pair<double[], double[]> getRule(BaseRuleFactory<? extends Number> factory,
int numberOfPoints)
throws NotStrictlyPositiveException, DimensionMismatchException {
return factory.getRule(numberOfPoints);
}
/**
* Performs a change of variable so that the integration can be performed
* on an arbitrary interval {@code [a, b]}.
* It is assumed that the natural interval is {@code [-1, 1]}.
*
* @param rule Original points and weights.
* @param a Lower bound of the integration interval.
* @param b Lower bound of the integration interval.
* @return the points and weights adapted to the new interval.
*/
private static Pair<double[], double[]> transform(Pair<double[], double[]> rule,
double a,
double b) {
final double[] points = rule.getFirst();
final double[] weights = rule.getSecond();
// Scaling
final double scale = (b - a) / 2;
final double shift = a + scale;
for (int i = 0; i < points.length; i++) {
points[i] = points[i] * scale + shift;
weights[i] *= scale;
}
return new Pair<double[], double[]>(points, weights);
}
}
|
