|
1 | 1 | package org.matheclipse.core.numerics.series.dp;
|
2 | 2 |
|
3 | 3 | /**
|
4 |
| - * Implements a modified adaptive Aitken delta^2 process for estimating infinite |
5 |
| - * series, based on the method in [1]. |
| 4 | + * Implements a modified adaptive Aitken delta^2 process for estimating infinite series, based on |
| 5 | + * the method in [1]. |
6 | 6 | *
|
7 | 7 | * <p>
|
8 | 8 | * References:
|
9 | 9 | * <ul>
|
10 |
| - * <li>[1] Osada, Naoki. Acceleration methods for slowly convergent sequences |
11 |
| - * and their applications. Diss. PhD thesis, Nagoya University, 1993.</li> |
| 10 | + * <li>[1] Osada, Naoki. Acceleration methods for slowly convergent sequences and their |
| 11 | + * applications. Diss. PhD thesis, Nagoya University, 1993.</li> |
12 | 12 | * </ul>
|
13 | 13 | * </p>
|
14 | 14 | */
|
15 | 15 | public final class Aitken extends SeriesAlgorithm {
|
16 | 16 |
|
17 |
| - private int status; |
18 |
| - private double xx, dx, dd, xp, alpha; |
19 |
| - private final double[] x; |
20 |
| - private final double[][] s, ds, ts, dts; |
| 17 | + private int status; |
| 18 | + private double xx, dx, dd, xp, alpha; |
| 19 | + private final double[] x; |
| 20 | + private final double[][] s, ds, ts, dts; |
21 | 21 |
|
22 |
| - public Aitken(final double tolerance, final int maxIters, final int patience) { |
23 |
| - super(tolerance, maxIters, patience); |
24 |
| - x = new double[maxIters]; |
25 |
| - s = new double[2][maxIters + 1]; |
26 |
| - ds = new double[2][maxIters + 1]; |
27 |
| - ts = new double[2][maxIters + 1]; |
28 |
| - dts = new double[2][maxIters + 1]; |
29 |
| - } |
| 22 | + public Aitken(final double tolerance, final int maxIters, final int patience) { |
| 23 | + super(tolerance, maxIters, patience); |
| 24 | + x = new double[maxIters]; |
| 25 | + s = new double[2][maxIters + 1]; |
| 26 | + ds = new double[2][maxIters + 1]; |
| 27 | + ts = new double[2][maxIters + 1]; |
| 28 | + dts = new double[2][maxIters + 1]; |
| 29 | + } |
30 | 30 |
|
31 |
| - @Override |
32 |
| - public final double next(final double e, final double term) { |
33 |
| - if (myIndex == 0) { |
34 |
| - xx = dx = dd = 0.0; |
35 |
| - } |
36 |
| - final double xxp = xx; |
37 |
| - x[myIndex] = xx = term; |
38 |
| - if (myIndex == 0) { |
39 |
| - dx = xx; |
40 |
| - } else if (myIndex == 1) { |
41 |
| - final double dxp = dx; |
42 |
| - dx = xx - xxp; |
43 |
| - dd = dx / (dx - dxp); |
44 |
| - } else { |
45 |
| - final double dxp = dx; |
46 |
| - final double ddp = dd; |
47 |
| - dx = xx - xxp; |
48 |
| - dd = dx / (dx - dxp); |
49 |
| - alpha = 1.0 / (dd - ddp) + 1.0; |
50 |
| - alpha = aitken(alpha, ts, dts, myIndex - 1, -2.0); |
51 |
| - } |
52 |
| - xp = xx; |
53 |
| - if (myIndex >= 2) { |
54 |
| - for (int m = 1; m <= myIndex + 1; ++m) { |
55 |
| - xp = aitken(x[m - 1], s, ds, m, alpha); |
56 |
| - } |
57 |
| - } |
58 |
| - ++myIndex; |
59 |
| - if (status == 1) { |
60 |
| - return Double.NaN; |
61 |
| - } else { |
62 |
| - return xp; |
63 |
| - } |
| 31 | + @Override |
| 32 | + public final double next(final double e, final double term) { |
| 33 | + if (myIndex == 0) { |
| 34 | + xx = dx = dd = 0.0; |
64 | 35 | }
|
65 |
| - |
66 |
| - private final double aitken(final double xx, final double[][] s, final double[][] ds, final int n, |
67 |
| - final double theta) { |
68 |
| - status = 0; |
69 |
| - if (n != 1) { |
70 |
| - final int kend = (n - 1) >> 1; |
71 |
| - System.arraycopy(s[1], 0, s[0], 0, kend + 1); |
72 |
| - if ((n & 1) == 0) { |
73 |
| - System.arraycopy(ds[1], 0, ds[0], 0, kend); |
74 |
| - } else { |
75 |
| - System.arraycopy(ds[1], 0, ds[0], 0, kend + 1); |
76 |
| - } |
77 |
| - } |
78 |
| - s[1][0] = xx; |
79 |
| - if (n == 1) { |
80 |
| - ds[1][0] = xx; |
81 |
| - return xx; |
82 |
| - } |
83 |
| - ds[1][0] = xx - s[0][0]; |
84 |
| - for (int k = 1; k <= (n >> 1); ++k) { |
85 |
| - final double w1 = ds[0][k - 1] * ds[1][k - 1]; |
86 |
| - final double w2 = ds[1][k - 1] - ds[0][k - 1]; |
87 |
| - if (Math.abs(w2) < TINY) { |
88 |
| - status = 1; |
89 |
| - return xx; |
90 |
| - } |
91 |
| - final int twok = k << 1; |
92 |
| - final double coef = ((twok - 1) - theta) / (twok - 2 - theta); |
93 |
| - s[1][k] = s[0][k - 1] - coef * w1 / w2; |
94 |
| - if (n != twok - 1) { |
95 |
| - ds[1][k] = s[1][k] - s[0][k]; |
96 |
| - } |
97 |
| - } |
98 |
| - final int kopt = n >> 1; |
99 |
| - return s[1][kopt]; |
| 36 | + final double xxp = xx; |
| 37 | + x[myIndex] = xx = term; |
| 38 | + if (myIndex == 0) { |
| 39 | + dx = xx; |
| 40 | + } else if (myIndex == 1) { |
| 41 | + final double dxp = dx; |
| 42 | + dx = xx - xxp; |
| 43 | + dd = dx / (dx - dxp); |
| 44 | + } else { |
| 45 | + final double dxp = dx; |
| 46 | + final double ddp = dd; |
| 47 | + dx = xx - xxp; |
| 48 | + dd = dx / (dx - dxp); |
| 49 | + alpha = 1.0 / (dd - ddp) + 1.0; |
| 50 | + alpha = aitken(alpha, ts, dts, myIndex - 1, -2.0); |
| 51 | + } |
| 52 | + xp = xx; |
| 53 | + if (myIndex >= 2) { |
| 54 | + for (int m = 1; m <= myIndex + 1; ++m) { |
| 55 | + xp = aitken(x[m - 1], s, ds, m, alpha); |
| 56 | + } |
| 57 | + } |
| 58 | + ++myIndex; |
| 59 | + if (status == 1) { |
| 60 | + return Double.NaN; |
| 61 | + } else { |
| 62 | + return xp; |
100 | 63 | }
|
| 64 | + } |
101 | 65 |
|
102 |
| - @Override |
103 |
| - public final String getName() { |
104 |
| - return "Modified Aitken"; |
| 66 | + private final double aitken(final double xx, final double[][] s, final double[][] ds, final int n, |
| 67 | + final double theta) { |
| 68 | + status = 0; |
| 69 | + if (n != 1) { |
| 70 | + final int kend = (n - 1) >> 1; |
| 71 | + System.arraycopy(s[1], 0, s[0], 0, kend + 1); |
| 72 | + if ((n & 1) == 0) { |
| 73 | + System.arraycopy(ds[1], 0, ds[0], 0, kend); |
| 74 | + } else { |
| 75 | + System.arraycopy(ds[1], 0, ds[0], 0, kend + 1); |
| 76 | + } |
105 | 77 | }
|
| 78 | + s[1][0] = xx; |
| 79 | + if (n == 1) { |
| 80 | + ds[1][0] = xx; |
| 81 | + return xx; |
| 82 | + } |
| 83 | + ds[1][0] = xx - s[0][0]; |
| 84 | + for (int k = 1; k <= (n >> 1); ++k) { |
| 85 | + final double w1 = ds[0][k - 1] * ds[1][k - 1]; |
| 86 | + final double w2 = ds[1][k - 1] - ds[0][k - 1]; |
| 87 | + if (Math.abs(w2) < TINY) { |
| 88 | + status = 1; |
| 89 | + return xx; |
| 90 | + } |
| 91 | + final int twok = k << 1; |
| 92 | + final double coef = ((twok - 1) - theta) / (twok - 2 - theta); |
| 93 | + s[1][k] = s[0][k - 1] - coef * w1 / w2; |
| 94 | + if (n != twok - 1) { |
| 95 | + ds[1][k] = s[1][k] - s[0][k]; |
| 96 | + } |
| 97 | + } |
| 98 | + final int kopt = n >> 1; |
| 99 | + return s[1][kopt]; |
| 100 | + } |
| 101 | + |
| 102 | + @Override |
| 103 | + public final String getName() { |
| 104 | + return "Modified Aitken"; |
| 105 | + } |
106 | 106 | }
|
0 commit comments