Skip to content

Commit eb6acdd

Browse files
axkraxkr
committed
WIP #1034 initial commit with modified library
- https://github.com/mike-gimelfarb/numerical-integration
1 parent 700c796 commit eb6acdd

25 files changed

+7112
-0
lines changed

‎symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/numerics/integral/ClenshawCurtis.java

+588
Large diffs are not rendered by default.

‎symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/numerics/integral/GaussKronrod.java

+1,195
Large diffs are not rendered by default.

‎symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/numerics/integral/GaussLegendre.java

+200
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,200 @@
1+
package org.matheclipse.core.numerics.integral;
2+
3+
import java.util.function.DoubleUnaryOperator;
4+
import org.matheclipse.core.numerics.utils.Constants;
5+
import org.matheclipse.core.numerics.utils.SimpleMath;
6+
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.
11+
*/
12+
public final class GaussLegendre extends Quadrature {
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);
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,
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);
199+
}
200+
}

‎symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/numerics/integral/GaussLobatto.java

+138
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,138 @@
1+
package org.matheclipse.core.numerics.integral;
2+
3+
import java.util.function.DoubleUnaryOperator;
4+
import org.matheclipse.core.numerics.utils.Constants;
5+
6+
/**
7+
* Implements an adaptive numerical integrator based on the 7-point
8+
* Gauss-Lobatto rule, as described in [1].
9+
*
10+
* <p>
11+
* References:
12+
* <ul>
13+
* <li>[1] Gander, W., Gautschi, W. Adaptive Quadrature�Revisited. BIT Numerical
14+
* Mathematics 40, 84�101 (2000). https://doi.org/10.1023/A:1022318402393</li>
15+
* </ul>
16+
* </p>
17+
*/
18+
public final class GaussLobatto extends Quadrature {
19+
20+
private static final double ALPHA = Constants.SQRT2 / Constants.SQRT3;
21+
private static final double BETA = 1.0 / Constants.SQRT5;
22+
23+
private static final double[] X = { 0.94288241569547971906, 0.64185334234578130578, 0.23638319966214988028 };
24+
private static final double[] Y = { 0.0158271919734801831, 0.0942738402188500455, 0.1550719873365853963,
25+
0.1888215739601824544, 0.1997734052268585268, 0.2249264653333395270, 0.2426110719014077338 };
26+
private static final double[] C = { 77.0, 432.0, 625.0, 672.0 };
27+
28+
private int fev;
29+
30+
public GaussLobatto(final double tolerance, final int maxEvaluations) {
31+
super(tolerance, maxEvaluations);
32+
}
33+
34+
@Override
35+
final QuadratureResult properIntegral(final DoubleUnaryOperator f, final double a,
36+
final double b) {
37+
return dlob8e(f, a, b);
38+
}
39+
40+
@Override
41+
public final String getName() {
42+
return "Gauss-Lobatto";
43+
}
44+
45+
private final QuadratureResult dlob8e(final DoubleUnaryOperator f, final double a,
46+
final double b) {
47+
48+
// compute interpolation points
49+
final double mid = 0.5 * (a + b);
50+
final double h = 0.5 * (b - a);
51+
final double y1 = f.applyAsDouble(a);
52+
final double y3 = f.applyAsDouble(mid - h * ALPHA);
53+
final double y5 = f.applyAsDouble(mid - h * BETA);
54+
final double y7 = f.applyAsDouble(mid);
55+
final double y9 = f.applyAsDouble(mid + h * BETA);
56+
final double y11 = f.applyAsDouble(mid + h * ALPHA);
57+
final double y13 = f.applyAsDouble(b);
58+
final double f1 = f.applyAsDouble(mid - h * X[0]);
59+
final double f2 = f.applyAsDouble(mid + h * X[0]);
60+
final double f3 = f.applyAsDouble(mid - h * X[1]);
61+
final double f4 = f.applyAsDouble(mid + h * X[1]);
62+
final double f5 = f.applyAsDouble(mid - h * X[2]);
63+
final double f6 = f.applyAsDouble(mid + h * X[2]);
64+
fev = 13;
65+
66+
// compute a crude initial estimate of the integral
67+
final double est1 = (y1 + y13 + 5.0 * (y5 + y9)) * (h / 6.0);
68+
69+
// compute a more refined estimate of the integral
70+
double est2 = C[0] * (y1 + y13) + C[1] * (y3 + y11) + C[2] * (y5 + y9) + C[3] * y7;
71+
est2 *= (h / 1470.0);
72+
73+
// compute the error estimate
74+
double s = Y[0] * (y1 + y13) + Y[1] * (f1 + f2);
75+
s += Y[2] * (y3 + y11) + Y[3] * (f3 + f4);
76+
s += Y[4] * (y5 + y9) + Y[5] * (f5 + f6) + Y[6] * y7;
77+
s *= h;
78+
double rtol = myTol;
79+
if (est1 != s) {
80+
final double r = Math.abs(est2 - s) / Math.abs(est1 - s);
81+
if (r > 0.0 && r < 1.0) {
82+
rtol /= r;
83+
}
84+
}
85+
double sign = Math.signum(s);
86+
if (sign == 0) {
87+
sign = 1.0;
88+
}
89+
double s1 = sign * Math.abs(s) * rtol / Constants.EPSILON;
90+
if (s == 0) {
91+
s1 = Math.abs(b - a);
92+
}
93+
94+
// call the recursive subroutine
95+
final double result = dlob8(f, a, b, y1, y13, s1, rtol);
96+
return new QuadratureResult(result, Double.NaN, fev, Double.isFinite(result));
97+
}
98+
99+
private final double dlob8(final DoubleUnaryOperator f, final double a, final double b,
100+
final double fa, final double fb, final double s, final double rtol) {
101+
102+
// check the budget of evaluations
103+
if (fev >= myMaxEvals) {
104+
return Double.NaN;
105+
}
106+
107+
// compute the interpolation points
108+
final double h = 0.5 * (b - a);
109+
final double mid = 0.5 * (a + b);
110+
final double mll = mid - h * ALPHA;
111+
final double ml = mid - h * BETA;
112+
final double mr = mid + BETA * h;
113+
final double mrr = mid + h * ALPHA;
114+
final double fmll = f.applyAsDouble(mll);
115+
final double fml = f.applyAsDouble(ml);
116+
final double fmid = f.applyAsDouble(mid);
117+
final double fmr = f.applyAsDouble(mr);
118+
final double fmrr = f.applyAsDouble(mrr);
119+
fev += 8;
120+
121+
// compute a crude estimate of the integral
122+
final double est1 = (fa + fb + 5.0 * (fml + fmr)) * (h / 6.0);
123+
124+
// compute a more refined estimate of the integral
125+
double est2 = C[0] * (fa + fb) + C[1] * (fmll + fmrr) + C[2] * (fml + fmr) + C[3] * fmid;
126+
est2 *= (h / 1470.0);
127+
128+
// check the convergence
129+
if (s + (est2 - est1) == s || mll <= a || b <= mrr) {
130+
return est2;
131+
}
132+
133+
// subdivide the integration region and repeat
134+
return dlob8(f, a, mll, fa, fmll, s, rtol) + dlob8(f, mll, ml, fmll, fml, s, rtol)
135+
+ dlob8(f, ml, mid, fml, fmid, s, rtol) + dlob8(f, mid, mr, fmid, fmr, s, rtol)
136+
+ dlob8(f, mr, mrr, fmr, fmrr, s, rtol) + dlob8(f, mrr, b, fmrr, fb, s, rtol);
137+
}
138+
}

0 commit comments

Comments
 (0)