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;
17
18 import com.irurueta.algebra.AlgebraException;
19 import com.irurueta.algebra.LUDecomposer;
20 import com.irurueta.algebra.Matrix;
21
22 /**
23 * Estimates the Padé approximant rational function by using a number of coefficients
24 * of a Taylor series.
25 * Padé approximants can yield more accurate solutions than Taylor series in certain situations.
26 */
27 public class PadeApproximantEstimator {
28
29 /**
30 * Number of times to iteratively improve LU decomposition solution by default.
31 */
32 private static final int DEFAULT_IMPROVEMENT_TIMES = 4;
33
34 /**
35 * Number of times to iteratively improve LU decomposition.
36 */
37 private final int improveTimes;
38
39 /**
40 * Computes LU decomposition to find denominator coefficients.
41 */
42 private final LUDecomposer luDecomposer = new LUDecomposer();
43
44 /**
45 * Number of coefficients being processed based on provided Taylor power series ones.
46 */
47 private int n;
48
49 /**
50 * Contains matrix to solve Padé coefficients.
51 */
52 private Matrix q;
53
54 /**
55 * Contains intermediate solution for denominator coefficients.
56 */
57 private Matrix x;
58
59 /**
60 * Contains intermediate solution for numerator coefficients.
61 */
62 private Matrix y;
63
64 /**
65 * Contains product of "a" and "x" matrices to iteratively improve LU solution.
66 */
67 private Matrix ax;
68
69 /**
70 * Contains residual to iteratively improve LU solution.
71 */
72 private Matrix residual;
73
74 /**
75 * Improved LU solution in one iteration.
76 */
77 private Matrix improvedX;
78
79 /**
80 * Default constructor.
81 * Uses default number of times to iteratively improve solution.
82 */
83 public PadeApproximantEstimator() {
84 this(DEFAULT_IMPROVEMENT_TIMES);
85 }
86
87 /**
88 * Constructor.
89 * @param improveTimes number of times to iteratively improve solution.
90 * @throws IllegalArgumentException if provided number of times is negative.
91 */
92 public PadeApproximantEstimator(final int improveTimes) {
93 if (improveTimes < 0) {
94 throw new IllegalArgumentException("Times must be zero or greater");
95 }
96
97 this.improveTimes = improveTimes;
98 }
99
100 /**
101 * Estimates Padé coefficients for provided Taylor power series ones.
102 *
103 * @param taylorCoefficients Taylor series coefficients.
104 * @return Result containing Padé approximant numerator and denominator coefficients.
105 * @throws NumericalException if a numerical error occurs.
106 * @throws IllegalArgumentException if provided number of Taylor series coefficients is less
107 * than 3.
108 */
109 public Result estimatePadeCoefficients(final double[] taylorCoefficients) throws NumericalException {
110 final var coefN = (taylorCoefficients.length - 1) / 2;
111 final var num = new double[coefN + 1];
112 final var denom = new double[coefN + 1];
113 estimatePadeCoefficients(taylorCoefficients, coefN, num, denom);
114 return new Result(num, denom);
115 }
116
117 /**
118 * Estimates Padé coefficients for provided Taylor power series ones.
119 *
120 * @param taylorCoefficients Taylor series coefficients.
121 * @param numeratorResult numerator coefficients of Padé approximant
122 * (must be (taylorCoefficients.length - 1) / 2).
123 * @param denominatorResult denominator coefficients of Padé approximant
124 * (must be (taylorCoefficients.length - 1) / 2)..
125 * @throws NumericalException if a numerical error occurs.
126 * @throws IllegalArgumentException if provided number of Taylor series coefficients is less
127 * than 3 or if provided numerator or denominator result coefficients have an invalid size.
128 */
129 public void estimatePadeCoefficients(
130 final double[] taylorCoefficients, final double[] numeratorResult, final double[] denominatorResult)
131 throws NumericalException {
132 final var coefN = (taylorCoefficients.length - 1) / 2;
133 estimatePadeCoefficients(taylorCoefficients, coefN, numeratorResult, denominatorResult);
134 }
135
136 /**
137 * Estimates Padé coefficients for provided Taylor power series ones.
138 *
139 * @param taylorCoefficients Taylor series coefficients.
140 * @param n Number of padé coefficients to generate.
141 * @param numeratorResult numerator coefficients of Padé approximant
142 * (must be (taylorCoefficients.length - 1) / 2).
143 * @param denominatorResult denominator coefficients of Padé approximant
144 * (must be (taylorCoefficients.length - 1) / 2)..
145 * @throws NumericalException if a numerical error occurs.
146 * @throws IllegalArgumentException if provided number of Taylor series coefficients is less
147 * than 3 or if provided numerator or denominator result coefficients have an invalid size.
148 */
149 private void estimatePadeCoefficients(
150 final double[] taylorCoefficients, final int n, final double[] numeratorResult,
151 final double[] denominatorResult) throws NumericalException {
152 if (taylorCoefficients.length < 3) {
153 throw new IllegalArgumentException("Length of Taylor series coefficients must be at least 3");
154 }
155
156 try {
157 // Based on Numerical Recipes section 5.12 Padé Approximants page 245.
158 final var nPlusOne = n +1;
159 if (numeratorResult.length != nPlusOne || denominatorResult.length != nPlusOne) {
160 throw new IllegalArgumentException("Wrong numerator or denominator array length");
161 }
162
163 if (this.n != n) {
164 initialize(n);
165 }
166 int j;
167 int k;
168 double sum;
169
170 for (j = 0; j < n; j++) {
171 // set up matrix for solving
172 y.setElementAtIndex(j, taylorCoefficients[n + j + 1]);
173 for (k = 0; k < n; k++) {
174 q.setElementAt(j, k, taylorCoefficients[j - k + n]);
175 }
176 }
177
178 luDecomposer.setInputMatrix(q);
179 luDecomposer.decompose();
180 luDecomposer.solve(y, x);
181
182 for (j = 0; j < improveTimes; j++) {
183 improveLuSolve(q, y, x, improvedX);
184 x.copyFrom(improvedX);
185 }
186
187 for (k = 0; k < n; k++) {
188 for (sum = taylorCoefficients[k + 1], j = 0; j <= k; j++) {
189 sum -= x.getElementAtIndex(j) * taylorCoefficients[k - j];
190 }
191 y.setElementAtIndex(k, sum);
192 }
193 numeratorResult[0] = taylorCoefficients[0];
194 denominatorResult[0] = 1.0;
195 for (j = 0; j < n; j++) {
196 numeratorResult[j + 1] = y.getElementAtIndex(j);
197 denominatorResult[j + 1] = -x.getElementAtIndex(j);
198 }
199
200 } catch (final AlgebraException ex) {
201 throw new NumericalException(ex);
202 }
203 }
204
205 /**
206 * One step to iteratively improve LU solve solution.
207 *
208 * @param a a matrix of a linear system of equations to be solved.
209 * @param b b matrix of a linear system of equations to be solved.
210 * @param x x matrix containing initial solution of linear system of equations to be improved.
211 * @param result matrix where result will be stored.
212 * @throws AlgebraException if a numerical error occurs.
213 */
214 private void improveLuSolve(final Matrix a, final Matrix b, final Matrix x, final Matrix result)
215 throws AlgebraException {
216 // Based on Numerical Recipes page 62
217 // We need to solve iteratively: A * deltaX = A * (x + deltaX) - b
218 // deltaX is the residual error between initially estimated x and the true x
219 // Hence:
220 // result = x - deltaX
221 // Where result will be closer to the true x
222
223 a.multiply(x, ax);
224 ax.subtract(b, residual);
225
226 luDecomposer.solve(residual, result);
227
228 result.multiplyByScalar(-1.0);
229 result.add(x);
230 }
231
232 /**
233 * Initializes required matrices.
234 *
235 * @param n length of required number of Padé coefficients for provided Taylor series ones.
236 * @throws AlgebraException if a numerical error occurs.
237 */
238 private void initialize(final int n) throws AlgebraException {
239 q = new Matrix(n, n);
240 x = new Matrix(n, 1);
241 y = new Matrix(n, 1);
242
243 if (improveTimes > 0) {
244 ax = new Matrix(n, 1);
245 residual = new Matrix(n, 1);
246 improvedX = new Matrix(n, 1);
247 }
248
249 this.n = n;
250 }
251
252 /**
253 * Contains result of Padé approximant.
254 */
255 public static class Result {
256
257 /**
258 * Numerator coefficients.
259 */
260 private final double[] numerators;
261
262 /**
263 * Denominator coefficients.
264 */
265 private final double[] denominators;
266
267 /**
268 * Constructor.
269 *
270 * @param numerators numerator coefficients.
271 * @param denominators denominator coefficients.
272 */
273 public Result(final double[] numerators, final double[] denominators) {
274 this.numerators = numerators;
275 this.denominators = denominators;
276 }
277
278 /**
279 * Gets numerator coefficients.
280 *
281 * @return numerator coefficients.
282 */
283 public double[] getNumerators() {
284 return numerators;
285 }
286
287 /**
288 * Gets denominator coefficients.
289 *
290 * @return denominator coefficients.
291 */
292 public double[] getDenominators() {
293 return denominators;
294 }
295 }
296 }