|
5 | 5 | import org.matheclipse.core.numerics.utils.SimpleMath;
|
6 | 6 |
|
7 | 7 | /**
|
8 |
| - * Integrate a real function of one variable over a finite interval using an |
9 |
| - * adaptive 8-point Legendre-Gauss algorithm. This algorithm is a translation of |
10 |
| - * the corresponding subroutine dgaus8 from the SLATEC library. |
| 8 | + * Integrate a real function of one variable over a finite interval using an adaptive 8-point |
| 9 | + * Legendre-Gauss algorithm. This algorithm is a translation of the corresponding subroutine dgaus8 |
| 10 | + * from the SLATEC library. |
11 | 11 | */
|
12 | 12 | public final class GaussLegendre extends Quadrature {
|
13 | 13 |
|
14 |
| - private static final double X1 = 1.83434642495649805e-1; |
15 |
| - private static final double X2 = 5.25532409916328986e-1; |
16 |
| - private static final double X3 = 7.96666477413626740e-1; |
17 |
| - private static final double X4 = 9.60289856497536232e-1; |
18 |
| - private static final double W1 = 3.62683783378361983e-1; |
19 |
| - private static final double W2 = 3.13706645877887287e-1; |
20 |
| - private static final double W3 = 2.22381034453374471e-1; |
21 |
| - private static final double W4 = 1.01228536290376259e-1; |
22 |
| - private static final double SQ2 = Constants.SQRT2; |
23 |
| - |
24 |
| - public GaussLegendre(final double tolerance, final int maxEvaluations) { |
25 |
| - super(tolerance, maxEvaluations); |
| 14 | + private static final double X1 = 1.83434642495649805e-1; |
| 15 | + private static final double X2 = 5.25532409916328986e-1; |
| 16 | + private static final double X3 = 7.96666477413626740e-1; |
| 17 | + private static final double X4 = 9.60289856497536232e-1; |
| 18 | + private static final double W1 = 3.62683783378361983e-1; |
| 19 | + private static final double W2 = 3.13706645877887287e-1; |
| 20 | + private static final double W3 = 2.22381034453374471e-1; |
| 21 | + private static final double W4 = 1.01228536290376259e-1; |
| 22 | + private static final double SQ2 = Constants.SQRT2; |
| 23 | + |
| 24 | + public GaussLegendre(final double tolerance, final int maxEvaluations) { |
| 25 | + super(tolerance, maxEvaluations); |
| 26 | + } |
| 27 | + |
| 28 | + @Override |
| 29 | + final QuadratureResult properIntegral(final DoubleUnaryOperator f, final double a, |
| 30 | + final double b) { |
| 31 | + |
| 32 | + // prepare variables |
| 33 | + final double[] err = {myTol}; |
| 34 | + final double[] ans = new double[1]; |
| 35 | + final int[] ierr = new int[1]; |
| 36 | + final int[] fev = new int[1]; |
| 37 | + |
| 38 | + // call main subroutine |
| 39 | + dgaus8(f, a, b, err, ans, ierr, fev); |
| 40 | + return new QuadratureResult(ans[0], err[0], fev[0], ierr[0] == 1 || ierr[0] == -1); |
| 41 | + } |
| 42 | + |
| 43 | + @Override |
| 44 | + public final String getName() { |
| 45 | + return "Gauss-Legendre"; |
| 46 | + } |
| 47 | + |
| 48 | + private void dgaus8(final DoubleUnaryOperator fun, final double a, final double b, |
| 49 | + final double[] err, final double[] ans, final int[] ierr, final int[] fev) { |
| 50 | + int k, kml = 6, kmx = 5000, l, lmx, mxl, nbits, nib, nlmx; |
| 51 | + double ae, anib, area, c, ce, ee, ef, eps, est, gl, glr, tol, vr; |
| 52 | + final int[] lr = new int[61]; |
| 53 | + final double[] aa = new double[61], hh = new double[61], vl = new double[61], |
| 54 | + gr = new double[61]; |
| 55 | + |
| 56 | + // Initialize |
| 57 | + fev[0] = 0; |
| 58 | + k = 53; |
| 59 | + anib = SimpleMath.D1MACH[5 - 1] * k / 0.30102000; |
| 60 | + nbits = (int) anib; |
| 61 | + nlmx = Math.min(60, (nbits * 5) / 8); |
| 62 | + ans[0] = 0.0; |
| 63 | + ierr[0] = 1; |
| 64 | + ce = 0.0; |
| 65 | + if (a == b) { |
| 66 | + if (err[0] < 0.0) { |
| 67 | + err[0] = ce; |
| 68 | + } |
| 69 | + return; |
26 | 70 | }
|
27 |
| - |
28 |
| - @Override |
29 |
| - final QuadratureResult properIntegral(final DoubleUnaryOperator f, final double a, |
30 |
| - final double b) { |
31 |
| - |
32 |
| - // prepare variables |
33 |
| - final double[] err = { myTol }; |
34 |
| - final double[] ans = new double[1]; |
35 |
| - final int[] ierr = new int[1]; |
36 |
| - final int[] fev = new int[1]; |
37 |
| - |
38 |
| - // call main subroutine |
39 |
| - dgaus8(f, a, b, err, ans, ierr, fev); |
40 |
| - return new QuadratureResult(ans[0], err[0], fev[0], ierr[0] == 1 || ierr[0] == -1); |
| 71 | + lmx = nlmx; |
| 72 | + if (b != 0.0 && SimpleMath.sign(1.0, b) * a > 0.0) { |
| 73 | + c = Math.abs(1.0 - a / b); |
| 74 | + if (c <= 0.1) { |
| 75 | + if (c <= 0.0) { |
| 76 | + if (err[0] < 0.0) { |
| 77 | + err[0] = ce; |
| 78 | + } |
| 79 | + return; |
| 80 | + } |
| 81 | + anib = 0.5 - Math.log(c) * Constants.LOG2_INV; |
| 82 | + nib = (int) anib; |
| 83 | + lmx = Math.min(nlmx, nbits - nib - 7); |
| 84 | + if (lmx < 1) { |
| 85 | + ierr[0] = -1; |
| 86 | + if (err[0] < 0.0) { |
| 87 | + err[0] = ce; |
| 88 | + } |
| 89 | + return; |
| 90 | + } |
| 91 | + } |
41 | 92 | }
|
42 | 93 |
|
43 |
| - @Override |
44 |
| - public final String getName() { |
45 |
| - return "Gauss-Legendre"; |
| 94 | + tol = Math.max(Math.abs(err[0]), SimpleMath.pow(2.0, 5 - nbits)) / 2.0; |
| 95 | + if (err[0] == 0.0) { |
| 96 | + tol = Math.sqrt(SimpleMath.D1MACH[4 - 1]); |
46 | 97 | }
|
47 |
| - |
48 |
| - private void dgaus8(final DoubleUnaryOperator fun, final double a, final double b, |
49 |
| - final double[] err, |
50 |
| - final double[] ans, final int[] ierr, final int[] fev) { |
51 |
| - int k, kml = 6, kmx = 5000, l, lmx, mxl, nbits, nib, nlmx; |
52 |
| - double ae, anib, area, c, ce, ee, ef, eps, est, gl, glr, tol, vr; |
53 |
| - final int[] lr = new int[61]; |
54 |
| - final double[] aa = new double[61], hh = new double[61], vl = new double[61], gr = new double[61]; |
55 |
| - |
56 |
| - // Initialize |
57 |
| - fev[0] = 0; |
58 |
| - k = 53; |
59 |
| - anib = SimpleMath.D1MACH[5 - 1] * k / 0.30102000; |
60 |
| - nbits = (int) anib; |
61 |
| - nlmx = Math.min(60, (nbits * 5) / 8); |
62 |
| - ans[0] = 0.0; |
63 |
| - ierr[0] = 1; |
64 |
| - ce = 0.0; |
65 |
| - if (a == b) { |
66 |
| - if (err[0] < 0.0) { |
67 |
| - err[0] = ce; |
68 |
| - } |
69 |
| - return; |
70 |
| - } |
71 |
| - lmx = nlmx; |
72 |
| - if (b != 0.0 && SimpleMath.sign(1.0, b) * a > 0.0) { |
73 |
| - c = Math.abs(1.0 - a / b); |
74 |
| - if (c <= 0.1) { |
75 |
| - if (c <= 0.0) { |
76 |
| - if (err[0] < 0.0) { |
77 |
| - err[0] = ce; |
78 |
| - } |
79 |
| - return; |
80 |
| - } |
81 |
| - anib = 0.5 - Math.log(c) * Constants.LOG2_INV; |
82 |
| - nib = (int) anib; |
83 |
| - lmx = Math.min(nlmx, nbits - nib - 7); |
84 |
| - if (lmx < 1) { |
85 |
| - ierr[0] = -1; |
86 |
| - if (err[0] < 0.0) { |
87 |
| - err[0] = ce; |
88 |
| - } |
89 |
| - return; |
90 |
| - } |
91 |
| - } |
92 |
| - } |
93 |
| - |
94 |
| - tol = Math.max(Math.abs(err[0]), SimpleMath.pow(2.0, 5 - nbits)) / 2.0; |
95 |
| - if (err[0] == 0.0) { |
96 |
| - tol = Math.sqrt(SimpleMath.D1MACH[4 - 1]); |
97 |
| - } |
98 |
| - eps = tol; |
99 |
| - hh[1 - 1] = (b - a) / 4.0; |
100 |
| - aa[1 - 1] = a; |
101 |
| - lr[1 - 1] = 1; |
102 |
| - l = 1; |
103 |
| - est = g8(fun, aa[l - 1] + 2.0 * hh[l - 1], 2.0 * hh[l - 1]); |
104 |
| - fev[0] += 8; |
105 |
| - k = 8; |
106 |
| - area = Math.abs(est); |
107 |
| - ef = 0.5; |
108 |
| - mxl = 0; |
109 |
| - |
110 |
| - while (true) { |
111 |
| - |
112 |
| - // Compute refined estimates, estimate the error, etc. |
113 |
| - if (fev[0] - 8 >= myMaxEvals) { |
114 |
| - // ans[0] = Double.NaN; |
115 |
| - ierr[0] = 2; |
116 |
| - return; |
117 |
| - } |
118 |
| - gl = g8(fun, aa[l - 1] + hh[l - 1], hh[l - 1]); |
119 |
| - fev[0] += 8; |
120 |
| - if (fev[0] - 8 >= myMaxEvals) { |
121 |
| - // ans[0] = Double.NaN; |
122 |
| - ierr[0] = 2; |
123 |
| - return; |
124 |
| - } |
125 |
| - gr[l - 1] = g8(fun, aa[l - 1] + 3.0 * hh[l - 1], hh[l - 1]); |
126 |
| - fev[0] += 8; |
127 |
| - |
128 |
| - k += 16; |
129 |
| - area += (Math.abs(gl) + Math.abs(gr[l - 1]) - Math.abs(est)); |
130 |
| - glr = gl + gr[l - 1]; |
131 |
| - ee = Math.abs(est - glr) * ef; |
132 |
| - ae = Math.max(eps * area, tol * Math.abs(glr)); |
133 |
| - if (ee - ae > 0.0) { |
134 |
| - |
135 |
| - // Consider the left half of this level |
136 |
| - if (k > kmx) { |
137 |
| - lmx = kml; |
138 |
| - } |
139 |
| - if (l >= lmx) { |
140 |
| - mxl = 1; |
141 |
| - } else { |
142 |
| - ++l; |
143 |
| - eps *= 0.5; |
144 |
| - ef /= SQ2; |
145 |
| - hh[l - 1] = hh[l - 1 - 1] * 0.5; |
146 |
| - lr[l - 1] = -1; |
147 |
| - aa[l - 1] = aa[l - 1 - 1]; |
148 |
| - est = gl; |
149 |
| - continue; |
150 |
| - } |
151 |
| - } |
152 |
| - |
153 |
| - ce += (est - glr); |
154 |
| - if (lr[l - 1] <= 0.0) { |
155 |
| - |
156 |
| - // Proceed to right half at this level |
157 |
| - vl[l - 1] = glr; |
158 |
| - est = gr[l - 1 - 1]; |
159 |
| - lr[l - 1] = 1; |
160 |
| - aa[l - 1] += 4.0 * hh[l - 1]; |
161 |
| - } else { |
162 |
| - |
163 |
| - // Return one level |
164 |
| - vr = glr; |
165 |
| - while (true) { |
166 |
| - if (l <= 1) { |
167 |
| - ans[0] = vr; |
168 |
| - if (mxl != 0 && Math.abs(ce) > 2.0 * tol * area) { |
169 |
| - ierr[0] = 2; |
170 |
| - } |
171 |
| - if (err[0] < 0.0) { |
172 |
| - err[0] = ce; |
173 |
| - } |
174 |
| - return; |
175 |
| - } |
176 |
| - --l; |
177 |
| - eps *= 2.0; |
178 |
| - ef *= SQ2; |
179 |
| - if (lr[l - 1] <= 0.0) { |
180 |
| - vl[l - 1] = vl[l + 1 - 1] + vr; |
181 |
| - est = gr[l - 1 - 1]; |
182 |
| - lr[l - 1] = 1; |
183 |
| - aa[l - 1] += 4.0 * hh[l - 1]; |
184 |
| - break; |
185 |
| - } else { |
186 |
| - vr += vl[l + 1 - 1]; |
187 |
| - } |
188 |
| - } |
189 |
| - } |
190 |
| - } |
191 |
| - } |
192 |
| - |
193 |
| - private static double g8(final DoubleUnaryOperator fun, final double x, final double h) { |
194 |
| - final double fx1 = fun.applyAsDouble(x - X1 * h) + fun.applyAsDouble(x + X1 * h); |
195 |
| - final double fx2 = fun.applyAsDouble(x - X2 * h) + fun.applyAsDouble(x + X2 * h); |
196 |
| - final double fx3 = fun.applyAsDouble(x - X3 * h) + fun.applyAsDouble(x + X3 * h); |
197 |
| - final double fx4 = fun.applyAsDouble(x - X4 * h) + fun.applyAsDouble(x + X4 * h); |
198 |
| - return h * (W1 * fx1 + W2 * fx2 + W3 * fx3 + W4 * fx4); |
| 98 | + eps = tol; |
| 99 | + hh[1 - 1] = (b - a) / 4.0; |
| 100 | + aa[1 - 1] = a; |
| 101 | + lr[1 - 1] = 1; |
| 102 | + l = 1; |
| 103 | + est = g8(fun, aa[l - 1] + 2.0 * hh[l - 1], 2.0 * hh[l - 1]); |
| 104 | + fev[0] += 8; |
| 105 | + k = 8; |
| 106 | + area = Math.abs(est); |
| 107 | + ef = 0.5; |
| 108 | + mxl = 0; |
| 109 | + |
| 110 | + while (true) { |
| 111 | + |
| 112 | + // Compute refined estimates, estimate the error, etc. |
| 113 | + if (fev[0] - 8 >= myMaxEvals) { |
| 114 | + // ans[0] = Double.NaN; |
| 115 | + ierr[0] = 2; |
| 116 | + return; |
| 117 | + } |
| 118 | + gl = g8(fun, aa[l - 1] + hh[l - 1], hh[l - 1]); |
| 119 | + fev[0] += 8; |
| 120 | + if (fev[0] - 8 >= myMaxEvals) { |
| 121 | + // ans[0] = Double.NaN; |
| 122 | + ierr[0] = 2; |
| 123 | + return; |
| 124 | + } |
| 125 | + gr[l - 1] = g8(fun, aa[l - 1] + 3.0 * hh[l - 1], hh[l - 1]); |
| 126 | + fev[0] += 8; |
| 127 | + |
| 128 | + k += 16; |
| 129 | + area += (Math.abs(gl) + Math.abs(gr[l - 1]) - Math.abs(est)); |
| 130 | + glr = gl + gr[l - 1]; |
| 131 | + ee = Math.abs(est - glr) * ef; |
| 132 | + ae = Math.max(eps * area, tol * Math.abs(glr)); |
| 133 | + if (ee - ae > 0.0) { |
| 134 | + |
| 135 | + // Consider the left half of this level |
| 136 | + if (k > kmx) { |
| 137 | + lmx = kml; |
| 138 | + } |
| 139 | + if (l >= lmx) { |
| 140 | + mxl = 1; |
| 141 | + } else { |
| 142 | + ++l; |
| 143 | + eps *= 0.5; |
| 144 | + ef /= SQ2; |
| 145 | + hh[l - 1] = hh[l - 1 - 1] * 0.5; |
| 146 | + lr[l - 1] = -1; |
| 147 | + aa[l - 1] = aa[l - 1 - 1]; |
| 148 | + est = gl; |
| 149 | + continue; |
| 150 | + } |
| 151 | + } |
| 152 | + |
| 153 | + ce += (est - glr); |
| 154 | + if (lr[l - 1] <= 0.0) { |
| 155 | + |
| 156 | + // Proceed to right half at this level |
| 157 | + vl[l - 1] = glr; |
| 158 | + est = gr[l - 1 - 1]; |
| 159 | + lr[l - 1] = 1; |
| 160 | + aa[l - 1] += 4.0 * hh[l - 1]; |
| 161 | + } else { |
| 162 | + |
| 163 | + // Return one level |
| 164 | + vr = glr; |
| 165 | + while (true) { |
| 166 | + if (l <= 1) { |
| 167 | + ans[0] = vr; |
| 168 | + if (mxl != 0 && Math.abs(ce) > 2.0 * tol * area) { |
| 169 | + ierr[0] = 2; |
| 170 | + } |
| 171 | + if (err[0] < 0.0) { |
| 172 | + err[0] = ce; |
| 173 | + } |
| 174 | + return; |
| 175 | + } |
| 176 | + --l; |
| 177 | + eps *= 2.0; |
| 178 | + ef *= SQ2; |
| 179 | + if (lr[l - 1] <= 0.0) { |
| 180 | + vl[l - 1] = vl[l + 1 - 1] + vr; |
| 181 | + est = gr[l - 1 - 1]; |
| 182 | + lr[l - 1] = 1; |
| 183 | + aa[l - 1] += 4.0 * hh[l - 1]; |
| 184 | + break; |
| 185 | + } else { |
| 186 | + vr += vl[l + 1 - 1]; |
| 187 | + } |
| 188 | + } |
| 189 | + } |
199 | 190 | }
|
| 191 | + } |
| 192 | + |
| 193 | + private static double g8(final DoubleUnaryOperator fun, final double x, final double h) { |
| 194 | + final double fx1 = fun.applyAsDouble(x - X1 * h) + fun.applyAsDouble(x + X1 * h); |
| 195 | + final double fx2 = fun.applyAsDouble(x - X2 * h) + fun.applyAsDouble(x + X2 * h); |
| 196 | + final double fx3 = fun.applyAsDouble(x - X3 * h) + fun.applyAsDouble(x + X3 * h); |
| 197 | + final double fx4 = fun.applyAsDouble(x - X4 * h) + fun.applyAsDouble(x + X4 * h); |
| 198 | + return h * (W1 * fx1 + W2 * fx2 + W3 * fx3 + W4 * fx4); |
| 199 | + } |
200 | 200 | }
|
0 commit comments