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 * Computes matrix function integration by using Romberg integration.
26 */
27 public class RombergTrapezoidalQuadratureMatrixIntegrator extends RombergMatrixIntegrator<TrapezoidalMatrixQuadrature> {
28
29 /**
30 * Default accuracy.
31 */
32 public static final double EPS = 1e-10;
33
34 /**
35 * Maximum number of allowed steps.
36 */
37 private static final int JMAX = 20;
38
39 /**
40 * Maximum number of allowed steps + 1.
41 */
42 private static final int JMAXP = JMAX + 1;
43
44 /**
45 * Minimum required number of steps.
46 */
47 private static final int K = 5;
48
49 /**
50 * Successive trapezoidal approximations.
51 */
52 private final Matrix[] s = new Matrix[JMAX];
53
54 /**
55 * Successive trapezoidal step sizes.
56 */
57 private final double[] h = new double[JMAXP];
58
59 /**
60 * Constructor.
61 *
62 * @param a Lower limit of integration.
63 * @param b Upper limit of integration.
64 * @param listener listener to evaluate a single dimension matrix (multivariate) function at
65 * required points.
66 * @param eps required accuracy.
67 * @throws WrongSizeException if size notified by provided listener is invalid.
68 */
69 public RombergTrapezoidalQuadratureMatrixIntegrator(
70 final double a, final double b, final MatrixSingleDimensionFunctionEvaluatorListener listener,
71 final double eps) throws WrongSizeException {
72 super(new TrapezoidalMatrixQuadrature(a, b, listener), eps);
73 }
74
75 /**
76 * Constructor with default accuracy.
77 *
78 * @param a Lower limit of integration.
79 * @param b Upper limit of integration.
80 * @param listener listener to evaluate a single dimension matrix (multivariate) function at
81 * required points.
82 * @throws WrongSizeException if size notified by provided listener is invalid.
83 */
84 public RombergTrapezoidalQuadratureMatrixIntegrator(
85 final double a, final double b, final MatrixSingleDimensionFunctionEvaluatorListener listener)
86 throws WrongSizeException {
87 this(a, b, listener, EPS);
88 }
89
90 /**
91 * Integrates function between provided lower and upper limits.
92 *
93 * @param result instance where result of integration will be stored.
94 * @throws IntegrationException if integration fails for numerical reasons.
95 */
96 @SuppressWarnings("Duplicates")
97 @Override
98 public void integrate(final Matrix result) throws IntegrationException {
99 try {
100 final var rows = q.getRows();
101 final var columns = q.getColumns();
102 final var elems = rows * columns;
103 for (var i = 0; i < JMAX; i++) {
104 s[i] = new Matrix(rows, columns);
105 }
106
107 final var interpolators = new PolynomialInterpolator[elems];
108 final var sInterp = new double[elems][JMAX];
109 for (var i = 0; i < elems; i++) {
110 sInterp[i] = new double[JMAX];
111 interpolators[i] = new PolynomialInterpolator(h, sInterp[i], K, false);
112 }
113
114 h[0] = 1.0;
115 for (var j = 1; j <= JMAX; j++) {
116 q.next(s[j - 1]);
117 // update sInterp
118 for (var i = 0; i < elems; i++) {
119 sInterp[i][j - 1] = s[j - 1].getElementAtIndex(i);
120 }
121 if (j >= K) {
122 var finished = true;
123 for (var i = 0; i < elems; i++) {
124 final var ss = interpolators[i].rawinterp(j - K, 0.0);
125 if (Double.isNaN(ss)) {
126 throw new IntegrationException("NaN was found");
127 }
128 result.setElementAtIndex(i, ss);
129 if (Math.abs(interpolators[i].getDy()) > eps * Math.abs(ss)) {
130 finished = false;
131 }
132 }
133
134 if (finished) {
135 return;
136 }
137 }
138 h[j] = 0.25 * h[j - 1];
139 }
140 } catch (final EvaluationException | InterpolationException | WrongSizeException e) {
141 throw new IntegrationException(e);
142 }
143
144 // Too many steps
145 throw new IntegrationException();
146 }
147
148 /**
149 * Gets type of quadrature.
150 *
151 * @return type of quadrature.
152 */
153 @Override
154 public QuadratureType getQuadratureType() {
155 return QuadratureType.TRAPEZOIDAL;
156 }
157 }