1 /*
2 * Copyright (C) 2023 Alberto Irurueta Carro (alberto@irurueta.com)
3 *
4 * Licensed under the Apache License, Version 2.0 (the "License");
5 * you may not use this file except in compliance with the License.
6 * You may obtain a copy of the License at
7 *
8 * http://www.apache.org/licenses/LICENSE-2.0
9 *
10 * Unless required by applicable law or agreed to in writing, software
11 * distributed under the License is distributed on an "AS IS" BASIS,
12 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13 * See the License for the specific language governing permissions and
14 * limitations under the License.
15 */
16 package com.irurueta.numerical.integration;
17
18 import com.irurueta.algebra.Matrix;
19 import com.irurueta.algebra.WrongSizeException;
20 import com.irurueta.numerical.EvaluationException;
21 import com.irurueta.numerical.interpolation.InterpolationException;
22 import com.irurueta.numerical.interpolation.PolynomialInterpolator;
23
24 /**
25 * Base integrator for implementations based on Romberg's method.
26 * Romberg's method is a generalization of Simpson's method for higher order integration schemes.
27 * This can be used to computed integration with less function evaluations for the same level of
28 * accuracy when more assumptions of function "smoothness" can be made.
29 * Implementations of Romberg's method are quite powerful for sufficiently smooth (e.g., analytic)
30 * integrands, integrated over intervals that contain no singularities, and where the endpoints are
31 * also non-singular. In such circumstances, Romberg's method, takes many, many fewer function
32 * evaluations than other method's such as Simpson's.
33 *
34 * @param <T> instance of a quadrature to be used for Romberg's method integration.
35 */
36 public abstract class RombergMatrixIntegrator<T extends MatrixQuadrature> extends MatrixIntegrator {
37
38 /**
39 * Default accuracy.
40 */
41 public static final double EPS = 3.0e-9;
42
43 /**
44 * Maximum number of allowed steps.
45 */
46 private static final int JMAX = 14;
47
48 /**
49 * Maximum number of allowed steps + 1.
50 */
51 private static final int JMAXP = JMAX + 1;
52
53 /**
54 * Minimum required number of steps.
55 */
56 private static final int K = 5;
57
58 /**
59 * Quadrature used for integration.
60 */
61 protected final T q;
62
63 /**
64 * Required accuracy.
65 */
66 protected final double eps;
67
68 /**
69 * Successive trapezoidal approximations.
70 */
71 private final Matrix[] s = new Matrix[JMAX];
72
73 /**
74 * Successive trapezoidal step sizes.
75 */
76 private final double[] h = new double[JMAXP];
77
78 /**
79 * Constructor.
80 *
81 * @param q Quadrature used for integration.
82 * @param eps Required accuracy.
83 */
84 protected RombergMatrixIntegrator(final T q, final double eps) {
85 this.q = q;
86 this.eps = eps;
87 }
88
89 /**
90 * Integrates function between provided lower and upper limits.
91 *
92 * @param result instance where result of integration will be stored.
93 * @throws IntegrationException if integration fails for numerical reasons.
94 */
95 @SuppressWarnings("Duplicates")
96 @Override
97 public void integrate(final Matrix result) throws IntegrationException {
98 try {
99 final var rows = q.getRows();
100 final var columns = q.getColumns();
101 final var elems = rows * columns;
102 for (var i = 0; i < JMAX; i++) {
103 s[i] = new Matrix(rows, columns);
104 }
105
106 final var interpolators = new PolynomialInterpolator[elems];
107 final var sInterp = new double[elems][JMAX];
108 for (int i = 0; i < elems; i++) {
109 sInterp[i] = new double[JMAX];
110 interpolators[i] = new PolynomialInterpolator(h, sInterp[i], K, false);
111 }
112
113 h[0] = 1.0;
114 for (int j = 1; j <= JMAX; j++) {
115 q.next(s[j - 1]);
116 // update sInterp
117 for (var i = 0; i < elems; i++) {
118 sInterp[i][j - 1] = s[j - 1].getElementAtIndex(i);
119 }
120 if (j >= K) {
121 var finished = true;
122 for (var i = 0; i < elems; i++) {
123 final var ss = interpolators[i].rawinterp(j - K, 0.0);
124 if (Double.isNaN(ss)) {
125 throw new IntegrationException("NaN was found");
126 }
127 result.setElementAtIndex(i, ss);
128 if (Math.abs(interpolators[i].getDy()) > eps * Math.abs(ss)) {
129 finished = false;
130 }
131 }
132
133 if (finished) {
134 return;
135 }
136 }
137 h[j] = h[j - 1] / 9.0;
138 }
139 } catch (final EvaluationException | InterpolationException | WrongSizeException e) {
140 throw new IntegrationException(e);
141 }
142
143 // Too many steps
144 throw new IntegrationException();
145 }
146
147 /**
148 * Gets type of integrator.
149 *
150 * @return type of integrator.
151 */
152 @Override
153 public IntegratorType getIntegratorType() {
154 return IntegratorType.ROMBERG;
155 }
156
157 /**
158 * Creates an integrator using Romberg's method.
159 * It must be noticed that upper limit of integration is ignored when using exponential
160 * mid-point quadrature type.
161 *
162 * @param a Lower limit of integration.
163 * @param b Upper limit of integration.
164 * @param listener listener to evaluate a single dimension function at required points.
165 * @param eps required accuracy.
166 * @param quadratureType quadrature type.
167 * @return created integrator.
168 * @throws IllegalArgumentException if provided quadrature type is not supported.
169 * @throws WrongSizeException if size notified by provided listener is invalid.
170 */
171 public static RombergMatrixIntegrator<MatrixQuadrature> create(
172 final double a, final double b, final MatrixSingleDimensionFunctionEvaluatorListener listener,
173 final double eps, final QuadratureType quadratureType) throws WrongSizeException {
174 return switch (quadratureType) {
175 case MID_POINT -> cast(new RombergMidPointQuadratureMatrixIntegrator(a, b, listener, eps));
176 case INFINITY_MID_POINT -> cast(new RombergInfinityMidPointQuadratureMatrixIntegrator(a, b, listener, eps));
177 case LOWER_SQUARE_ROOT_MID_POINT ->
178 cast(new RombergLowerSquareRootMidPointQuadratureMatrixIntegrator(a, b, listener, eps));
179 case UPPER_SQUARE_ROOT_MID_POINT ->
180 cast(new RombergUpperSquareRootMidPointQuadratureMatrixIntegrator(a, b, listener, eps));
181 case EXPONENTIAL_MID_POINT ->
182 cast(new RombergExponentialMidPointQuadratureMatrixIntegrator(a, listener, eps));
183 case DOUBLE_EXPONENTIAL_RULE ->
184 cast(new RombergDoubleExponentialRuleQuadratureMatrixIntegrator(a, b, listener, eps));
185 default -> cast(new RombergTrapezoidalQuadratureMatrixIntegrator(a, b, listener, eps));
186 };
187 }
188
189 /**
190 * Creates an integrator using Romberg's method and having default accuracy.
191 * It must be noticed that upper limit of integration is ignored when using exponential
192 * mid-point quadrature type.
193 *
194 * @param a Lower limit of integration.
195 * @param b Upper limit of integration.
196 * @param listener listener to evaluate a single dimension function at required points.
197 * @param quadratureType quadrature type.
198 * @return created integrator.
199 * @throws IllegalArgumentException if provided quadrature type is not supported.
200 * @throws WrongSizeException if size notified by provided listener is invalid.
201 */
202 public static RombergMatrixIntegrator<MatrixQuadrature> create(
203 final double a, final double b, final MatrixSingleDimensionFunctionEvaluatorListener listener,
204 final QuadratureType quadratureType) throws WrongSizeException {
205 return switch (quadratureType) {
206 case MID_POINT -> cast(new RombergMidPointQuadratureMatrixIntegrator(a, b, listener));
207 case INFINITY_MID_POINT -> cast(new RombergInfinityMidPointQuadratureMatrixIntegrator(a, b, listener));
208 case LOWER_SQUARE_ROOT_MID_POINT ->
209 cast(new RombergLowerSquareRootMidPointQuadratureMatrixIntegrator(a, b, listener));
210 case UPPER_SQUARE_ROOT_MID_POINT ->
211 cast(new RombergUpperSquareRootMidPointQuadratureMatrixIntegrator(a, b, listener));
212 case EXPONENTIAL_MID_POINT -> cast(new RombergExponentialMidPointQuadratureMatrixIntegrator(a, listener));
213 case DOUBLE_EXPONENTIAL_RULE ->
214 cast(new RombergDoubleExponentialRuleQuadratureMatrixIntegrator(a, b, listener));
215 default -> cast(new RombergTrapezoidalQuadratureMatrixIntegrator(a, b, listener));
216 };
217 }
218
219 /**
220 * Creates an integrator using Romberg's method and default quadrature type.
221 * It must be noticed that upper limit of integration is ignored when using exponential
222 * mid-point quadrature type.
223 *
224 * @param a Lower limit of integration.
225 * @param b Upper limit of integration.
226 * @param listener listener to evaluate a single dimension function at required points.
227 * @param eps required accuracy.
228 * @return created integrator.
229 * @throws WrongSizeException if size notified by provided listener is invalid.
230 */
231 public static RombergMatrixIntegrator<MatrixQuadrature> create(
232 final double a, final double b, final MatrixSingleDimensionFunctionEvaluatorListener listener,
233 final double eps) throws WrongSizeException {
234 return create(a, b, listener, eps, DEFAULT_QUADRATURE_TYPE);
235 }
236
237 /**
238 * Creates an integrator using Romberg's method and having default accuracy and quadrature type.
239 * It must be noticed that upper limit of integration is ignored when using exponential
240 * mid-point quadrature type.
241 *
242 * @param a Lower limit of integration.
243 * @param b Upper limit of integration.
244 * @param listener listener to evaluate a single dimension function at required points.
245 * @return created integrator.
246 * @throws WrongSizeException if size notified by provided listener is invalid.
247 */
248 public static RombergMatrixIntegrator<MatrixQuadrature> create(
249 final double a, final double b, final MatrixSingleDimensionFunctionEvaluatorListener listener)
250 throws WrongSizeException {
251 return create(a, b, listener, DEFAULT_QUADRATURE_TYPE);
252 }
253
254 /**
255 * Casts integrator to a quadrature integrator without wildcard parameter.
256 *
257 * @param integrator integrator to be cast.
258 * @return cast integrator.
259 */
260 private static RombergMatrixIntegrator<MatrixQuadrature> cast(final RombergMatrixIntegrator<?> integrator) {
261 //noinspection unchecked
262 return (RombergMatrixIntegrator<MatrixQuadrature>) integrator;
263 }
264 }