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