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.numerical.EvaluationException;
19 import com.irurueta.numerical.SingleDimensionFunctionEvaluatorListener;
20 import com.irurueta.numerical.interpolation.InterpolationException;
21 import com.irurueta.numerical.interpolation.PolynomialInterpolator;
22
23 /**
24 * Base integrator for implementations based on Romberg's method.
25 * Romberg's method is a generalization of Simpson's method for higher order integration schemes.
26 * This can be used to computed integration with less function evaluations for the same level of
27 * accuracy when more assumptions of function "smoothness" can be made.
28 * Implementations of Romberg's method are quite powerful for sufficiently smooth (e.g., analytic)
29 * integrands, integrated over intervals that contain no singularities, and where the endpoints are
30 * also non-singular. In such circumstances, Romberg's method, takes many, many fewer function
31 * evaluations than other method's such as Simpson's.
32 *
33 * @param <T> instance of a quadrature to be used for Romber'gs method integration.
34 */
35 public abstract class RombergIntegrator<T extends Quadrature> extends Integrator {
36
37 /**
38 * Default accuracy.
39 */
40 public static final double EPS = 3.0e-9;
41
42 /**
43 * Maximum number of allowed steps.
44 */
45 private static final int JMAX = 14;
46
47 /**
48 * Maximum number of allowed steps + 1.
49 */
50 private static final int JMAXP = JMAX + 1;
51
52 /**
53 * Minimum required number of steps.
54 */
55 private static final int K = 5;
56
57 /**
58 * Quadrature used for integration.
59 */
60 protected final T q;
61
62 /**
63 * Required accuracy.
64 */
65 protected final double eps;
66
67 /**
68 * Successive trapezoidal approximations.
69 */
70 private final double[] s = new double[JMAX];
71
72 /**
73 * Successive trapezoidal step sizes.
74 */
75 private final double[] h = new double[JMAXP];
76
77 /**
78 * Polynomial interpolator.
79 */
80 private final PolynomialInterpolator interpolator = new PolynomialInterpolator(h, s, K, false);
81
82 /**
83 * Constructor.
84 *
85 * @param q Quadrature used for integration.
86 * @param eps Required accuracy.
87 */
88 protected RombergIntegrator(final T q, final double eps) {
89 this.q = q;
90 this.eps = eps;
91 }
92
93 /**
94 * Integrates function between lower and upper limits defined by provided quadrature.
95 * This implementation is suitable for open intervals.
96 *
97 * @return result of integration.
98 * @throws IntegrationException if integration fails for numerical reasons.
99 */
100 @Override
101 public double integrate() throws IntegrationException {
102 try {
103 h[0] = 1.0;
104 for (var j = 1; j <= JMAX; j++) {
105 s[j - 1] = q.next();
106 if (j >= K) {
107 final var ss = interpolator.rawinterp(j - K, 0.0);
108 if (Math.abs(interpolator.getDy()) <= eps * Math.abs(ss)) {
109 return ss;
110 }
111 }
112 h[j] = h[j - 1] / 9.0;
113 }
114 } catch (final EvaluationException | InterpolationException e) {
115 throw new IntegrationException(e);
116 }
117
118 // Too many steps
119 throw new IntegrationException();
120 }
121
122 /**
123 * Gets type of integrator.
124 *
125 * @return type of integrator.
126 */
127 @Override
128 public IntegratorType getIntegratorType() {
129 return IntegratorType.ROMBERG;
130 }
131
132 /**
133 * Creates an integrator using Romberg's method.
134 * It must be noticed that upper limit of integration is ignored when using exponential
135 * mid-point quadrature type.
136 *
137 * @param a Lower limit of integration.
138 * @param b Upper limit of integration.
139 * @param listener listener to evaluate a single dimension function at required points.
140 * @param eps required accuracy.
141 * @param quadratureType quadrature type.
142 * @return created integrator.
143 */
144 public static RombergIntegrator<Quadrature> create(
145 final double a, final double b, final SingleDimensionFunctionEvaluatorListener listener, final double eps,
146 final QuadratureType quadratureType) {
147 return switch (quadratureType) {
148 case MID_POINT -> cast(new RombergMidPointQuadratureIntegrator(a, b, listener, eps));
149 case INFINITY_MID_POINT -> cast(new RombergInfinityMidPointQuadratureIntegrator(a, b, listener, eps));
150 case LOWER_SQUARE_ROOT_MID_POINT ->
151 cast(new RombergLowerSquareRootMidPointQuadratureIntegrator(a, b, listener, eps));
152 case UPPER_SQUARE_ROOT_MID_POINT ->
153 cast(new RombergUpperSquareRootMidPointQuadratureIntegrator(a, b, listener, eps));
154 case EXPONENTIAL_MID_POINT -> cast(new RombergExponentialMidPointQuadratureIntegrator(a, listener, eps));
155 case DOUBLE_EXPONENTIAL_RULE ->
156 cast(new RombergDoubleExponentialRuleQuadratureIntegrator(a, b, listener, eps));
157 default -> cast(new RombergTrapezoidalQuadratureIntegrator(a, b, listener, eps));
158 };
159 }
160
161 /**
162 * Creates an integrator using Romberg's method and having default accuracy.
163 * It must be noticed that upper limit of integration is ignored when using exponential
164 * mid-point quadrature type.
165 *
166 * @param a Lower limit of integration.
167 * @param b Upper limit of integration.
168 * @param listener listener to evaluate a single dimension function at required points.
169 * @param quadratureType quadrature type.
170 * @return created integrator.
171 */
172 public static RombergIntegrator<Quadrature> create(
173 final double a, final double b, final SingleDimensionFunctionEvaluatorListener listener,
174 final QuadratureType quadratureType) {
175 return switch (quadratureType) {
176 case MID_POINT -> cast(new RombergMidPointQuadratureIntegrator(a, b, listener));
177 case INFINITY_MID_POINT -> cast(new RombergInfinityMidPointQuadratureIntegrator(a, b, listener));
178 case LOWER_SQUARE_ROOT_MID_POINT ->
179 cast(new RombergLowerSquareRootMidPointQuadratureIntegrator(a, b, listener));
180 case UPPER_SQUARE_ROOT_MID_POINT ->
181 cast(new RombergUpperSquareRootMidPointQuadratureIntegrator(a, b, listener));
182 case EXPONENTIAL_MID_POINT -> cast(new RombergExponentialMidPointQuadratureIntegrator(a, listener));
183 case DOUBLE_EXPONENTIAL_RULE -> cast(new RombergDoubleExponentialRuleQuadratureIntegrator(a, b, listener));
184 default -> cast(new RombergTrapezoidalQuadratureIntegrator(a, b, listener));
185 };
186 }
187
188 /**
189 * Creates an integrator using Romberg's method and default quadrature type.
190 * It must be noticed that upper limit of integration is ignored when using exponential
191 * mid-point quadrature type.
192 *
193 * @param a Lower limit of integration.
194 * @param b Upper limit of integration.
195 * @param listener listener to evaluate a single dimension function at required points.
196 * @param eps required accuracy.
197 * @return created integrator.
198 */
199 public static RombergIntegrator<Quadrature> create(
200 final double a, final double b, final SingleDimensionFunctionEvaluatorListener listener, final double eps) {
201 return create(a, b, listener, eps, DEFAULT_QUADRATURE_TYPE);
202 }
203
204 /**
205 * Creates an integrator using Romberg's method and having default accuracy and quadrature type.
206 * It must be noticed that upper limit of integration is ignored when using exponential
207 * mid-point quadrature type.
208 *
209 * @param a Lower limit of integration.
210 * @param b Upper limit of integration.
211 * @param listener listener to evaluate a single dimension function at required points.
212 * @return created integrator.
213 */
214 public static RombergIntegrator<Quadrature> create(
215 final double a, final double b, final SingleDimensionFunctionEvaluatorListener listener) {
216 return create(a, b, listener, DEFAULT_QUADRATURE_TYPE);
217 }
218
219 /**
220 * Casts integrator to a quadrature integrator without wildcard parameter.
221 * .
222 * @param integrator integrator to be cast.
223 * @return cast integrator.
224 */
225 private static RombergIntegrator<Quadrature> cast(final RombergIntegrator<?> integrator) {
226 //noinspection unchecked
227 return (RombergIntegrator<Quadrature>) integrator;
228 }
229 }