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 matrix quadrature using trapezoidal algorithm.
24 * This implementation is suitable for non-improper integrands, which consist
25 * of functions with not known singularities that can be evaluated on all the integration
26 * interval, which must be finite.
27 */
28 public class TrapezoidalMatrixQuadrature extends MatrixQuadrature {
29
30 /**
31 * Lower limit of integration.
32 */
33 private final double a;
34
35 /**
36 * Upper limit of integration.
37 */
38 private final double b;
39
40 /**
41 * Number of rows of quadrature result.
42 */
43 private final int rows;
44
45 /**
46 * Number of columns of quadrature result.
47 */
48 private final int columns;
49
50 /**
51 * Current value of integral.
52 */
53 private final Matrix s;
54
55 /**
56 * Listener to evaluate single dimension matrix functions at required points.
57 */
58 private final MatrixSingleDimensionFunctionEvaluatorListener listener;
59
60 /**
61 * Temporary value storing evaluation at point "a".
62 */
63 private final Matrix tmpA;
64
65 /**
66 * Temporary value storing evaluation at point "b".
67 */
68 private final Matrix tmpB;
69
70 /**
71 * Temporary value storing summation of evaluations.
72 */
73 private final Matrix sum;
74
75 /**
76 * Temporary value storing evaluation at point x.
77 */
78 private final Matrix tmpX;
79
80 /**
81 * Constructor.
82 *
83 * @param a Lower limit of integration.
84 * @param b Upper limit of integration.
85 * @param listener listener to evaluate a single dimension function at required points.
86 * @throws WrongSizeException if size notified by provided listener is invalid.
87 */
88 public TrapezoidalMatrixQuadrature(
89 final double a, final double b, final MatrixSingleDimensionFunctionEvaluatorListener listener)
90 throws WrongSizeException {
91 this.n = 0;
92 this.a = a;
93 this.b = b;
94
95 rows = listener.getRows();
96 columns = listener.getColumns();
97 s = new Matrix(rows, columns);
98 tmpA = new Matrix(rows, columns);
99 tmpB = new Matrix(rows, columns);
100 sum = new Matrix(rows, columns);
101 tmpX = new Matrix(rows, columns);
102
103 this.listener = listener;
104 }
105
106 /**
107 * Gets lower limit of integration.
108 *
109 * @return lower limit of integration.
110 */
111 public double getA() {
112 return a;
113 }
114
115 /**
116 * Gets upper limit of integration.
117 *
118 * @return upper limit of integration.
119 */
120 public double getB() {
121 return b;
122 }
123
124 /**
125 * Gets current value of integral.
126 *
127 * @return current vlaue of integral.
128 */
129 public Matrix getS() {
130 return s;
131 }
132
133 /**
134 * Returns the value of the integral at the nth stage of refinement.
135 *
136 * @param result instance where the value of the integral at the nth stage of refinement will
137 * be stored.
138 * @throws EvaluationException Raised if something failed during the evaluation.
139 */
140 @Override
141 public void next(final Matrix result) throws EvaluationException {
142 try {
143 double x;
144 double tnm;
145 double del;
146 int it;
147 int j;
148 n++;
149 if (n == 1) {
150 // (s = 0.5 * (b - a) * (listener.evaluate(a) + listener.evaluate(b))
151 listener.evaluate(a, tmpA);
152 listener.evaluate(b, tmpB);
153 s.copyFrom(tmpA);
154 s.add(tmpB);
155 s.multiplyByScalar(0.5 * (b - a));
156 result.copyFrom(s);
157 } else {
158 for (it = 1, j = 1; j < n - 1; j++) {
159 it <<= 1;
160 }
161 tnm = it;
162 // This is the spacing of the points to be added
163 del = (b - a) / tnm;
164 x = a + 0.5 * del;
165 sum.initialize(0.0);
166 for (j = 0; j < it; j++, x += del) {
167 listener.evaluate(x, tmpX);
168 sum.add(tmpX);
169 }
170 // This replaces s by its refined value
171 // s = 0.5 * (s + (b-a) * sum / tnm)
172 sum.multiplyByScalar((b - a) / tnm);
173 s.add(sum);
174 s.multiplyByScalar(0.5);
175 result.copyFrom(s);
176 }
177 } catch (final WrongSizeException ex) {
178 throw new EvaluationException(ex);
179 }
180 }
181
182 /**
183 * Gets type of quadrature.
184 *
185 * @return type of quadrature.
186 */
187 @Override
188 public QuadratureType getType() {
189 return QuadratureType.TRAPEZOIDAL;
190 }
191
192 /**
193 * Gets number of rows of quadrature result.
194 *
195 * @return number of rows of quadrature result.
196 */
197 @Override
198 protected int getRows() {
199 return rows;
200 }
201
202 /**
203 * Gets number of columns of quadrature result.
204 *
205 * @return number of columns of quadrature result.
206 */
207 @Override
208 protected int getColumns() {
209 return columns;
210 }
211 }