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
22 /**
23 * Implementation of quadrature using double exponential, which allows integration with a variable
24 * transformation.
25 * This implementation is suitable for improper integrands containing singularities.
26 */
27 public class DoubleExponentialRuleMatrixQuadrature extends MatrixQuadrature {
28
29 /**
30 * Default transformation of range of integration.
31 */
32 public static final double DEFAULT_HMAX = 3.7;
33
34 /**
35 * Lower limit of integration.
36 */
37 private final double a;
38
39 /**
40 * Upper limit of integration.
41 */
42 private final double b;
43
44 /**
45 * Maximum step size. Determines transformation of range of integration.
46 */
47 private final double hmax;
48
49 /**
50 * Number of rows of quadrature result.
51 */
52 private final int rows;
53
54 /**
55 * Number of columns of quadrature result.
56 */
57 private final int columns;
58
59 /**
60 * Value of the next stage of refinement.
61 */
62 private final Matrix s;
63
64 /**
65 * Temporary value storing summation of evaluations.
66 */
67 private final Matrix sum;
68
69 /**
70 * Temporary value storing evaluation at mid-point.
71 */
72 private final Matrix tmpMid;
73
74 /**
75 * Temporary value storing evaluation at lower bound.
76 */
77 private final Matrix tmpA;
78
79 /**
80 * Temporary value storing evaluation at upper bound.
81 */
82 private final Matrix tmpB;
83
84 /**
85 * Temporary value storing evaluation at point x.
86 */
87 private final Matrix tmpX;
88
89 /**
90 * Listener to evaluate single dimension functions at required points.
91 */
92 private final DoubleExponentialMatrixSingleDimensionFunctionEvaluatorListener listener;
93
94 /**
95 * Constructor.
96 *
97 * @param listener listener to evaluate a single dimension matrix function at required points.
98 * @param a Lower limit of integration.
99 * @param b Upper limit of integration.
100 * @param hmax Maximum step size. This quadrature transforms the range of integration to
101 * [-hmax, hmax].
102 * @throws WrongSizeException if size notified by provided listener is invalid.
103 */
104 public DoubleExponentialRuleMatrixQuadrature(
105 final DoubleExponentialMatrixSingleDimensionFunctionEvaluatorListener listener, final double a,
106 final double b, final double hmax) throws WrongSizeException {
107 this.listener = listener;
108 this.a = a;
109 this.b = b;
110 this.hmax = hmax;
111 n = 0;
112
113 rows = listener.getRows();
114 columns = listener.getColumns();
115 s = new Matrix(rows, columns);
116 sum = new Matrix(rows, columns);
117 tmpMid = new Matrix(rows, columns);
118 tmpA = new Matrix(rows, columns);
119 tmpB = new Matrix(rows, columns);
120 tmpX = new Matrix(rows, columns);
121 }
122
123 /**
124 * Constructor.
125 *
126 * @param listener listener to evaluate a single dimension matrix function at required points.
127 * @param a Lower limit of integration.
128 * @param b Upper limit of integration.
129 * @throws WrongSizeException if size notified by provided listener is invalid.
130 */
131 public DoubleExponentialRuleMatrixQuadrature(
132 final DoubleExponentialMatrixSingleDimensionFunctionEvaluatorListener listener, final double a,
133 final double b) throws WrongSizeException {
134 this(listener, a, b, DEFAULT_HMAX);
135 }
136
137 /**
138 * Constructor.
139 *
140 * @param listener listener to evaluate a single dimension matrix function at required points.
141 * @param a Lower limit of integration.
142 * @param b Upper limit of integration.
143 * @param hmax Maximum step size. This quadrature transforms the range of integration to
144 * [-hmax, hmax].
145 * @throws WrongSizeException if size notified by provided listener is invalid.
146 */
147 public DoubleExponentialRuleMatrixQuadrature(
148 final MatrixSingleDimensionFunctionEvaluatorListener listener, final double a, final double b,
149 final double hmax) throws WrongSizeException {
150 this(new DoubleExponentialMatrixSingleDimensionFunctionEvaluatorListener() {
151 @Override
152 public void evaluate(double x, double delta, Matrix result) throws EvaluationException {
153 listener.evaluate(x, result);
154 }
155
156 @Override
157 public int getRows() {
158 return listener.getRows();
159 }
160
161 @Override
162 public int getColumns() {
163 return listener.getColumns();
164 }
165 }, a, b, hmax);
166 }
167
168 /**
169 * Constructor.
170 *
171 * @param listener listener to evaluate a single dimension matrix function at required points.
172 * @param a Lower limit of integration.
173 * @param b Upper limit of integration.
174 * @throws WrongSizeException if size notified by provided listener is invalid.
175 */
176 public DoubleExponentialRuleMatrixQuadrature(
177 final MatrixSingleDimensionFunctionEvaluatorListener listener, final double a, final double b)
178 throws WrongSizeException {
179 this(listener, a, b, DEFAULT_HMAX);
180 }
181
182
183 /**
184 * Returns the value of the integral at the nth stage of refinement.
185 *
186 * @param result instance where the value of the integral at the nth stage of refinement will
187 * be stored.
188 * @throws EvaluationException Raised if something failed during the evaluation.
189 */
190 @SuppressWarnings("Duplicates")
191 @Override
192 public void next(Matrix result) throws EvaluationException {
193 try {
194 // On the first call to the function next (n = 1), the routine returns the crudest estimate
195 // of integral between a and b of f(x). Subsequent calls to next (n = 2, 3, ...) will
196 // improve the accuracy by adding 2^(n-1) additional interior points.
197 double del;
198 double fact;
199 double q;
200 double t;
201 final double twoh;
202 int it;
203 int j;
204 n++;
205 if (n == 1) {
206 fact = 0.25;
207 // s = hmax * 2.0 * (b - a) * fact * listener.evaluate(0.5 * (b + a), 0.5 * (b - a))
208 listener.evaluate(0.5 * (b + a), 0.5 * (b - a), tmpMid);
209 s.copyFrom(tmpMid);
210 s.multiplyByScalar(hmax * 2.0 * (b - a) * fact);
211 result.copyFrom(s);
212 } else {
213 for (it = 1, j = 1; j < n - 1; j++) {
214 it <<= 1;
215 }
216
217 // Twice the spacing of the points to be added
218 twoh = hmax / it;
219 t = 0.5 * twoh;
220 sum.initialize(0.0);
221 for (j = 0; j < it; j++) {
222 q = Math.exp(-2.0 * Math.sinh(t));
223 del = (b - a) * q / (1.0 + q);
224 final var value = 1.0 + q;
225 fact = q / (value * value) * Math.cosh(t);
226 listener.evaluate(a + del, del, tmpA);
227 listener.evaluate(b - del, del, tmpB);
228 tmpX.copyFrom(tmpA);
229 tmpX.add(tmpB);
230 tmpX.multiplyByScalar(fact);
231 sum.add(tmpX);
232 t += twoh;
233 }
234
235 // Replace s by its refined value and return.
236 // s = 0.5 * s + (b - a) * twoh * sum
237 sum.multiplyByScalar((b - a) * twoh);
238 s.multiplyByScalar(0.5);
239 s.add(sum);
240 result.copyFrom(s);
241 }
242 } catch (final WrongSizeException ex) {
243 throw new EvaluationException(ex);
244 }
245 }
246
247 /**
248 * Gets type of quadrature.
249 *
250 * @return type of quadrature.
251 */
252 @Override
253 public QuadratureType getType() {
254 return QuadratureType.DOUBLE_EXPONENTIAL_RULE;
255 }
256
257 /**
258 * Gets number of rows of quadrature result.
259 *
260 * @return number of rows of quadrature result.
261 */
262 @Override
263 protected int getRows() {
264 return rows;
265 }
266
267 /**
268 * Gets number of columns of quadrature result.
269 *
270 * @return number of columns of quadrature result.
271 */
272 @Override
273 protected int getColumns() {
274 return columns;
275 }
276 }