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.AlgebraException;
19 import com.irurueta.algebra.Matrix;
20 import com.irurueta.algebra.WrongSizeException;
21 import com.irurueta.numerical.EvaluationException;
22
23 /**
24 * Base integrator for implementations based on Simpson's method.
25 * Simpson's method is an optimization of Trapezoidal quadrature integrator.
26 * Implementations of this class will in general be more efficient than
27 * Trapezoidal quadrature matrix integrators (i.e., require fewer function evaluations) when the
28 * matrix function to be integrated has a finite fourth derivative (i.e., a continuous third
29 * derivative).
30 *
31 * @param <T> a quadrature.
32 */
33 public abstract class SimpsonMatrixIntegrator<T extends MatrixQuadrature> extends MatrixIntegrator {
34
35 /**
36 * Default accuracy.
37 */
38 public static final double EPS = 1e-10;
39
40 /**
41 * Minimum required number of steps.
42 */
43 private static final int JMIN = 5;
44
45 /**
46 * Maximum number of allowed steps.
47 */
48 private static final int JMAX = 20;
49
50 /**
51 * Quadrature used for integration.
52 */
53 private final T q;
54
55 /**
56 * Required accuracy.
57 */
58 private final double eps;
59
60 /**
61 * Constructor.
62 *
63 * @param q Quadrature used for integration.
64 * @param eps Required accuracy.
65 */
66 protected SimpsonMatrixIntegrator(final T q, final double eps) {
67 this.q = q;
68 this.eps = eps;
69 }
70
71 /**
72 * Integrates function between provided lower and upper limits.
73 *
74 * @param result instance where result of integration will be stored.
75 * @throws IntegrationException if integration fails for numerical reasons.
76 */
77 @SuppressWarnings("Duplicates")
78 @Override
79 public void integrate(final Matrix result) throws IntegrationException {
80 try {
81 final var rows = q.getRows();
82 final var columns = q.getColumns();
83 final var st = new Matrix(rows, columns);
84 final var s = new Matrix(rows, columns);
85 final var ost = new Matrix(rows, columns);
86 final var os = new Matrix(rows, columns);
87 final var tmp = new Matrix(rows, columns);
88 for (var j = 0; j < JMAX; j++) {
89 q.next(st);
90
91 // s = (4.0 * st - ost) / 3.0
92 tmp.copyFrom(st);
93 tmp.multiplyByScalar(4.0);
94 tmp.subtract(ost);
95 tmp.multiplyByScalar(1.0 / 3.0);
96 s.copyFrom(tmp);
97
98 if (j > JMIN && (Math.abs(normMin(s) - normMin(os)) < eps * normMin(os)
99 || (normMin(s) == 0.0 && normMin(os) == 0.0))) {
100 result.copyFrom(s);
101 return;
102 }
103 os.copyFrom(s);
104 ost.copyFrom(st);
105 }
106 } catch (final EvaluationException | AlgebraException e) {
107 throw new IntegrationException(e);
108 }
109
110 // Too many steps
111 throw new IntegrationException();
112 }
113
114 /**
115 * Gets type of integrator.
116 *
117 * @return type of integrator.
118 */
119 @Override
120 public IntegratorType getIntegratorType() {
121 return IntegratorType.SIMPSON;
122 }
123
124 /**
125 * Creates an integrator using Simpson's method.
126 *
127 * @param a Lower limit of integration.
128 * @param b Upper limit of integration.
129 * @param listener listener to evaluate a single dimension function at required points.
130 * @param eps required accuracy.
131 * @param quadratureType quadrature type.
132 * @return created integrator.
133 * @throws IllegalArgumentException if provided quadrature type is not supported.
134 * @throws WrongSizeException if size notified by provided listener is invalid.
135 */
136 public static SimpsonMatrixIntegrator<MatrixQuadrature> create(
137 final double a, final double b, final MatrixSingleDimensionFunctionEvaluatorListener listener,
138 final double eps, final QuadratureType quadratureType) throws WrongSizeException {
139 return switch (quadratureType) {
140 case TRAPEZOIDAL -> cast(new SimpsonTrapezoidalQuadratureMatrixIntegrator(a, b, listener, eps));
141 case MID_POINT -> cast(new SimpsonMidPointQuadratureMatrixIntegrator(a, b, listener, eps));
142 case INFINITY_MID_POINT -> cast(new SimpsonInfinityMidPointQuadratureMatrixIntegrator(a, b, listener, eps));
143 case LOWER_SQUARE_ROOT_MID_POINT ->
144 cast(new SimpsonLowerSquareRootMidPointQuadratureMatrixIntegrator(a, b, listener, eps));
145 case UPPER_SQUARE_ROOT_MID_POINT ->
146 cast(new SimpsonUpperSquareRootMidPointQuadratureMatrixIntegrator(a, b, listener, eps));
147 case DOUBLE_EXPONENTIAL_RULE ->
148 cast(new SimpsonDoubleExponentialRuleQuadratureMatrixIntegrator(a, b, listener, eps));
149 default -> throw new IllegalArgumentException();
150 };
151 }
152
153 /**
154 * Creates an integrator using Simpson's method and having default accuracy.
155 *
156 * @param a Lower limit of integration.
157 * @param b Upper limit of integration.
158 * @param listener listener to evaluate a single dimension function at required points.
159 * @param quadratureType quadrature type.
160 * @return created integrator.
161 * @throws IllegalArgumentException if provided quadrature type is not supported.
162 * @throws WrongSizeException if size notified by provided listener is invalid.
163 */
164 public static SimpsonMatrixIntegrator<MatrixQuadrature> create(
165 final double a, final double b, final MatrixSingleDimensionFunctionEvaluatorListener listener,
166 final QuadratureType quadratureType) throws WrongSizeException {
167 return switch (quadratureType) {
168 case TRAPEZOIDAL -> cast(new SimpsonTrapezoidalQuadratureMatrixIntegrator(a, b, listener));
169 case MID_POINT -> cast(new SimpsonMidPointQuadratureMatrixIntegrator(a, b, listener));
170 case INFINITY_MID_POINT -> cast(new SimpsonInfinityMidPointQuadratureMatrixIntegrator(a, b, listener));
171 case LOWER_SQUARE_ROOT_MID_POINT ->
172 cast(new SimpsonLowerSquareRootMidPointQuadratureMatrixIntegrator(a, b, listener));
173 case UPPER_SQUARE_ROOT_MID_POINT ->
174 cast(new SimpsonUpperSquareRootMidPointQuadratureMatrixIntegrator(a, b, listener));
175 case DOUBLE_EXPONENTIAL_RULE ->
176 cast(new SimpsonDoubleExponentialRuleQuadratureMatrixIntegrator(a, b, listener));
177 default -> throw new IllegalArgumentException();
178 };
179 }
180
181 /**
182 * Creates an integrator using Simpson's method and default quadrature type.
183 *
184 * @param a Lower limit of integration.
185 * @param b Upper limit of integration.
186 * @param listener listener to evaluate a single dimension function at required points.
187 * @param eps required accuracy.
188 * @return created integrator.
189 * @throws WrongSizeException if size notified by provided listener is invalid.
190 */
191 public static SimpsonMatrixIntegrator<MatrixQuadrature> create(
192 final double a, final double b, final MatrixSingleDimensionFunctionEvaluatorListener listener,
193 final double eps) throws WrongSizeException {
194 return create(a, b, listener, eps, DEFAULT_QUADRATURE_TYPE);
195 }
196
197 /**
198 * Creates an integrator using Simpson's method and having default accuracy and quadrature type.
199 *
200 * @param a Lower limit of integration.
201 * @param b Upper limit of integration.
202 * @param listener listener to evaluate a single dimension function at required points.
203 * @return created integrator.
204 * @throws WrongSizeException if size notified by provided listener is invalid.
205 */
206 public static SimpsonMatrixIntegrator<MatrixQuadrature> create(
207 final double a, final double b, final MatrixSingleDimensionFunctionEvaluatorListener listener)
208 throws WrongSizeException {
209 return create(a, b, listener, DEFAULT_QUADRATURE_TYPE);
210 }
211
212 /**
213 * Estimates smallest norm of provided matrix.
214 * Smallest norm is used to ensure convergence of all elements in matrix.
215 *
216 * @param a matrix to compute min norm for.
217 * @return estimated min norm.
218 */
219 private static double normMin(final Matrix a) {
220 var min = Double.MAX_VALUE;
221 var buffer = a.getBuffer();
222 for (var v : buffer) {
223 final var value = Math.abs(v);
224 if (Double.isNaN(value)) {
225 return value;
226 }
227
228 if (value < min) {
229 min = value;
230 }
231 }
232 return min;
233 }
234
235 /**
236 * Casts integrator to a quadrature integrator without wildcard parameter.
237 *
238 * @param integrator integrator to be cast.
239 * @return cast integrator.
240 */
241 private static SimpsonMatrixIntegrator<MatrixQuadrature> cast(final SimpsonMatrixIntegrator<?> integrator) {
242 //noinspection unchecked
243 return (SimpsonMatrixIntegrator<MatrixQuadrature>) integrator;
244 }
245 }