1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.apache.commons.math.analysis.integration;
18
19 import org.apache.commons.math.ConvergenceException;
20 import org.apache.commons.math.FunctionEvaluationException;
21 import org.apache.commons.math.MathRuntimeException;
22 import org.apache.commons.math.MaxIterationsExceededException;
23 import org.apache.commons.math.analysis.UnivariateRealFunction;
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 public class LegendreGaussIntegrator extends UnivariateRealIntegratorImpl {
53
54
55 private static final double[] ABSCISSAS_2 = {
56 -1.0 / Math.sqrt(3.0),
57 1.0 / Math.sqrt(3.0)
58 };
59
60
61 private static final double[] WEIGHTS_2 = {
62 1.0,
63 1.0
64 };
65
66
67 private static final double[] ABSCISSAS_3 = {
68 -Math.sqrt(0.6),
69 0.0,
70 Math.sqrt(0.6)
71 };
72
73
74 private static final double[] WEIGHTS_3 = {
75 5.0 / 9.0,
76 8.0 / 9.0,
77 5.0 / 9.0
78 };
79
80
81 private static final double[] ABSCISSAS_4 = {
82 -Math.sqrt((15.0 + 2.0 * Math.sqrt(30.0)) / 35.0),
83 -Math.sqrt((15.0 - 2.0 * Math.sqrt(30.0)) / 35.0),
84 Math.sqrt((15.0 - 2.0 * Math.sqrt(30.0)) / 35.0),
85 Math.sqrt((15.0 + 2.0 * Math.sqrt(30.0)) / 35.0)
86 };
87
88
89 private static final double[] WEIGHTS_4 = {
90 (90.0 - 5.0 * Math.sqrt(30.0)) / 180.0,
91 (90.0 + 5.0 * Math.sqrt(30.0)) / 180.0,
92 (90.0 + 5.0 * Math.sqrt(30.0)) / 180.0,
93 (90.0 - 5.0 * Math.sqrt(30.0)) / 180.0
94 };
95
96
97 private static final double[] ABSCISSAS_5 = {
98 -Math.sqrt((35.0 + 2.0 * Math.sqrt(70.0)) / 63.0),
99 -Math.sqrt((35.0 - 2.0 * Math.sqrt(70.0)) / 63.0),
100 0.0,
101 Math.sqrt((35.0 - 2.0 * Math.sqrt(70.0)) / 63.0),
102 Math.sqrt((35.0 + 2.0 * Math.sqrt(70.0)) / 63.0)
103 };
104
105
106 private static final double[] WEIGHTS_5 = {
107 (322.0 - 13.0 * Math.sqrt(70.0)) / 900.0,
108 (322.0 + 13.0 * Math.sqrt(70.0)) / 900.0,
109 128.0 / 225.0,
110 (322.0 + 13.0 * Math.sqrt(70.0)) / 900.0,
111 (322.0 - 13.0 * Math.sqrt(70.0)) / 900.0
112 };
113
114
115 private final double[] abscissas;
116
117
118 private final double[] weights;
119
120
121
122
123
124
125
126 public LegendreGaussIntegrator(final int n, final int defaultMaximalIterationCount)
127 throws IllegalArgumentException {
128 super(defaultMaximalIterationCount);
129 switch(n) {
130 case 2 :
131 abscissas = ABSCISSAS_2;
132 weights = WEIGHTS_2;
133 break;
134 case 3 :
135 abscissas = ABSCISSAS_3;
136 weights = WEIGHTS_3;
137 break;
138 case 4 :
139 abscissas = ABSCISSAS_4;
140 weights = WEIGHTS_4;
141 break;
142 case 5 :
143 abscissas = ABSCISSAS_5;
144 weights = WEIGHTS_5;
145 break;
146 default :
147 throw MathRuntimeException.createIllegalArgumentException(
148 "{0} points Legendre-Gauss integrator not supported, " +
149 "number of points must be in the {1}-{2} range",
150 n, 2, 5);
151 }
152
153 }
154
155
156 @Deprecated
157 public double integrate(final double min, final double max)
158 throws ConvergenceException, FunctionEvaluationException, IllegalArgumentException {
159 return integrate(f, min, max);
160 }
161
162
163 public double integrate(final UnivariateRealFunction f,
164 final double min, final double max)
165 throws ConvergenceException, FunctionEvaluationException, IllegalArgumentException {
166
167 clearResult();
168 verifyInterval(min, max);
169 verifyIterationCount();
170
171
172 double oldt = stage(f, min, max, 1);
173
174 int n = 2;
175 for (int i = 0; i < maximalIterationCount; ++i) {
176
177
178 final double t = stage(f, min, max, n);
179
180
181 final double delta = Math.abs(t - oldt);
182 final double limit =
183 Math.max(absoluteAccuracy,
184 relativeAccuracy * (Math.abs(oldt) + Math.abs(t)) * 0.5);
185
186
187 if ((i + 1 >= minimalIterationCount) && (delta <= limit)) {
188 setResult(t, i);
189 return result;
190 }
191
192
193 double ratio = Math.min(4, Math.pow(delta / limit, 0.5 / abscissas.length));
194 n = Math.max((int) (ratio * n), n + 1);
195 oldt = t;
196
197 }
198
199 throw new MaxIterationsExceededException(maximalIterationCount);
200
201 }
202
203
204
205
206
207
208
209
210
211
212
213 private double stage(final UnivariateRealFunction f,
214 final double min, final double max, final int n)
215 throws FunctionEvaluationException {
216
217
218 final double step = (max - min) / n;
219 final double halfStep = step / 2.0;
220
221
222 double midPoint = min + halfStep;
223 double sum = 0.0;
224 for (int i = 0; i < n; ++i) {
225 for (int j = 0; j < abscissas.length; ++j) {
226 sum += weights[j] * f.value(midPoint + halfStep * abscissas[j]);
227 }
228 midPoint += step;
229 }
230
231 return halfStep * sum;
232
233 }
234
235 }