This is an automated email from the ASF dual-hosted git repository. asf-gitbox-commits pushed a commit to branch geoapi-4.0 in repository https://gitbox.apache.org/repos/asf/sis.git
commit 1b20c15dd36e066e142a307de75fecd4409a2adb Author: jsorel <[email protected]> AuthorDate: Tue Oct 13 12:23:47 2026 +0200 feat(Geometry): add GaussLegendreQuadrature --- .../operation/GaussLegendreQuadrature.java | 196 +++++++++++++++++++++ .../operation/GaussLegendreQuadratureTest.java | 194 ++++++++++++++++++++ 2 files changed, 390 insertions(+) diff --git a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/operation/GaussLegendreQuadrature.java b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/operation/GaussLegendreQuadrature.java new file mode 100644 index 0000000000..fdeee31810 --- /dev/null +++ b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/operation/GaussLegendreQuadrature.java @@ -0,0 +1,196 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +package org.apache.sis.geometries.operation; + +import java.util.function.DoubleUnaryOperator; + + +/** + * Numerical integration of a function of one variable by the Gauss-Legendre quadrature. + * + * <p>A rule of <var>n</var> points evaluates the function at the roots of the Legendre polynomial of + * degree <var>n</var>, mapped to the integration interval, and is exact for polynomials of degree up + * to 2<var>n</var> − 1. For smooth functions the error decreases very fast with <var>n</var>. For + * functions which are not smooth everywhere, the {@linkplain #integrate(DoubleUnaryOperator, double, + * double, double, int, IntervalConsumer) adaptive integration} halves the interval where needed.</p> + * + * <p>Instances are immutable and thread-safe.</p> + * + * @author Johann Sorel (Geomatys) + * @see <a href="https://en.wikipedia.org/wiki/Gauss%E2%80%93Legendre_quadrature";>Gauss–Legendre quadrature</a> + */ +public final class GaussLegendreQuadrature { + + /** + * Receives the sub-intervals over which an adaptive integration has converged, + * in increasing order of the integration variable. + */ + @FunctionalInterface + public interface IntervalConsumer { + /** + * Receives a sub-interval and the integral of the function over it. + * + * @param start start of the sub-interval. + * @param end end of the sub-interval. + * @param integral integral of the function from {@code start} to {@code end}. + */ + void accept(double start, double end, double integral); + } + + /** + * Abscissas of the rule in the [−1 … +1] interval, in increasing order. + */ + private final double[] nodes; + + /** + * Weight of each abscissa. + */ + private final double[] weights; + + /** + * Creates a rule of the given number of points. + * + * @param n number of points, at least 1. + * @throws IllegalArgumentException if the number of points is less than 1. + */ + public GaussLegendreQuadrature(final int n) { + if (n < 1) { + throw new IllegalArgumentException("The number of points must be at least 1, but got " + n + '.'); + } + nodes = new double[n]; + weights = new double[n]; + /* + * The roots are symmetric around zero, so only the positive ones are computed, by Newton + * iterations starting from an asymptotic approximation of the root. The Legendre polynomial + * and its derivative are evaluated by the three-term recurrence. + */ + for (int i = 0; i < (n + 1) / 2; i++) { + double z = Math.cos(Math.PI * (i + 0.75) / (n + 0.5)); + double derivative; + for (int iteration = 0; ; iteration++) { + double p1 = 1, p2 = 0; + for (int j = 1; j <= n; j++) { + final double p3 = p2; + p2 = p1; + p1 = ((2*j - 1) * z * p2 - (j - 1) * p3) / j; + } + derivative = n * (z * p1 - p2) / (z * z - 1); + final double previous = z; + z -= p1 / derivative; + if (Math.abs(z - previous) <= Math.ulp(1.0) || iteration >= 100) break; + } + final double weight = 2 / ((1 - z * z) * derivative * derivative); + nodes[i] = -z; + nodes[n - 1 - i] = z; + weights[i] = weight; + weights[n - 1 - i] = weight; + } + if ((n & 1) != 0) { + nodes[n / 2] = 0; // Exact zero instead of a value close to it. + } + } + + /** + * Returns the number of points of this rule. + * + * @return number of points. + */ + public int getNumPoints() { + return nodes.length; + } + + /** + * Returns the abscissas of this rule in the [−1 … +1] interval, in increasing order. + * + * @return a copy of the abscissas. + */ + public double[] getNodes() { + return nodes.clone(); + } + + /** + * Returns the weights of this rule, in the order of the {@linkplain #getNodes() abscissas}. + * + * @return a copy of the weights. + */ + public double[] getWeights() { + return weights.clone(); + } + + /** + * Integrates the given function from {@code a} to {@code b} by a single application of this rule. + * The result is negative if {@code b} is less than {@code a} and the function is positive. + * + * @param f the function to integrate. + * @param a lower bound of the integral. + * @param b upper bound of the integral. + * @return integral of the function from {@code a} to {@code b}. + */ + public double integrate(final DoubleUnaryOperator f, final double a, final double b) { + final double half = (b - a) / 2; + final double mid = (a + b) / 2; + double sum = 0; + for (int i = 0; i < nodes.length; i++) { + sum += weights[i] * f.applyAsDouble(mid + half * nodes[i]); + } + return sum * half; + } + + /** + * Integrates the given function from {@code a} to {@code b}, halving the interval until the + * integrals over the two halves agree with the integral over the whole interval within the given + * relative tolerance. Each sub-interval is halved independently, so the subdivision is finer only + * where the function needs it. + * + * <p>If the function cannot be integrated within the tolerance, for example because its integral + * is zero or because it is not continuous, the halving stops after {@code maxDepth} levels.</p> + * + * @param f the function to integrate. + * @param a lower bound of the integral. + * @param b upper bound of the integral. + * @param tolerance relative tolerance on the integral of each sub-interval. + * @param maxDepth maximal number of interval halvings. + * @param consumer receives the sub-intervals and their integrals in order, or {@code null} if none. + * @return integral of the function from {@code a} to {@code b}. + */ + public double integrate(final DoubleUnaryOperator f, final double a, final double b, + final double tolerance, final int maxDepth, final IntervalConsumer consumer) + { + return adaptive(f, a, b, integrate(f, a, b), tolerance, maxDepth, consumer); + } + + /** + * Recursive part of the adaptive integration, given the integral over the whole interval. + */ + private double adaptive(final DoubleUnaryOperator f, final double a, final double b, final double whole, + final double tolerance, final int depth, final IntervalConsumer consumer) + { + final double m = (a + b) / 2; + final double left = integrate(f, a, m); + final double right = integrate(f, m, b); + final double sum = left + right; + if (depth <= 0 || Math.abs(sum - whole) <= tolerance * Math.abs(sum)) { + if (consumer != null) { + consumer.accept(a, m, left); + consumer.accept(m, b, right); + } + return sum; + } + return adaptive(f, a, m, left, tolerance, depth - 1, consumer) + + adaptive(f, m, b, right, tolerance, depth - 1, consumer); + } +} diff --git a/incubator/src/org.apache.sis.geometry/test/org/apache/sis/geometries/operation/GaussLegendreQuadratureTest.java b/incubator/src/org.apache.sis.geometry/test/org/apache/sis/geometries/operation/GaussLegendreQuadratureTest.java new file mode 100644 index 0000000000..f5fd2ab678 --- /dev/null +++ b/incubator/src/org.apache.sis.geometry/test/org/apache/sis/geometries/operation/GaussLegendreQuadratureTest.java @@ -0,0 +1,194 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +package org.apache.sis.geometries.operation; + +import java.util.ArrayList; +import java.util.List; + +// Test dependencies +import static org.junit.jupiter.api.Assertions.*; +import org.junit.jupiter.api.Test; + + +/** + * Tests {@link GaussLegendreQuadrature}. + * + * @author Johann Sorel (Geomatys) + */ +public class GaussLegendreQuadratureTest { + /** + * Tolerance on the computed abscissas and weights. + */ + private static final double EPS = 1E-15; + + @Test + public void testInvalidNumberOfPoints() { + assertThrows(IllegalArgumentException.class, () -> new GaussLegendreQuadrature(0)); + assertThrows(IllegalArgumentException.class, () -> new GaussLegendreQuadrature(-1)); + } + + /** + * The rules of 1 and 2 points have simple closed forms. + */ + @Test + public void testSmallRules() { + final GaussLegendreQuadrature one = new GaussLegendreQuadrature(1); + assertEquals(1, one.getNumPoints()); + assertArrayEquals(new double[] {0}, one.getNodes(), EPS); + assertArrayEquals(new double[] {2}, one.getWeights(), EPS); + + final GaussLegendreQuadrature two = new GaussLegendreQuadrature(2); + final double x = 1 / Math.sqrt(3); + assertArrayEquals(new double[] {-x, x}, two.getNodes(), EPS); + assertArrayEquals(new double[] { 1, 1}, two.getWeights(), EPS); + } + + /** + * Compares the rule of 5 points with its closed form. + */ + @Test + public void testFivePoints() { + final GaussLegendreQuadrature rule = new GaussLegendreQuadrature(5); + final double x1 = Math.sqrt(5 - 2 * Math.sqrt(10.0 / 7)) / 3; + final double x2 = Math.sqrt(5 + 2 * Math.sqrt(10.0 / 7)) / 3; + final double w1 = (322 + 13 * Math.sqrt(70)) / 900; + final double w2 = (322 - 13 * Math.sqrt(70)) / 900; + assertArrayEquals(new double[] {-x2, -x1, 0, x1, x2}, rule.getNodes(), EPS); + assertArrayEquals(new double[] {w2, w1, 128.0 / 225, w1, w2}, rule.getWeights(), EPS); + } + + /** + * Compares the rule of 10 points with the values tabulated in the literature. + */ + @Test + public void testTenPoints() { + final GaussLegendreQuadrature rule = new GaussLegendreQuadrature(10); + final double[] x = {0.1488743389816312, 0.4333953941292472, 0.6794095682990244, 0.8650633666889845, 0.9739065285171717}; + final double[] w = {0.2955242247147529, 0.2692667193099963, 0.2190863625159820, 0.1494513491505806, 0.0666713443086881}; + final double[] nodes = rule.getNodes(); + final double[] weights = rule.getWeights(); + for (int i = 0; i < x.length; i++) { + assertEquals( x[i], nodes [5 + i], EPS); + assertEquals(-x[i], nodes [4 - i], EPS); + assertEquals( w[i], weights[5 + i], EPS); + assertEquals( w[i], weights[4 - i], EPS); + } + } + + /** + * For any number of points, the abscissas are increasing and symmetric around zero, + * and the weights are positive and sum to 2, the length of the [−1 … +1] interval. + */ + @Test + public void testNodesAndWeightsProperties() { + for (int n = 1; n <= 64; n++) { + final GaussLegendreQuadrature rule = new GaussLegendreQuadrature(n); + final double[] nodes = rule.getNodes(); + final double[] weights = rule.getWeights(); + double sum = 0; + for (int i = 0; i < n; i++) { + if (i > 0) { + assertTrue(nodes[i] > nodes[i-1], "Abscissas shall be increasing."); + } + assertTrue(nodes[i] > -1 && nodes[i] < 1, "Abscissas shall be inside the interval."); + assertEquals(-nodes[i], nodes[n - 1 - i], EPS, "Abscissas shall be symmetric."); + assertTrue(weights[i] > 0, "Weights shall be positive."); + sum += weights[i]; + } + assertEquals(2, sum, 1E-13, "Weights shall sum to 2."); + } + } + + /** + * A rule of <var>n</var> points is exact for polynomials of degree up to 2<var>n</var> − 1, + * but not beyond. + */ + @Test + public void testPolynomialExactness() { + for (int n = 1; n <= 12; n++) { + final GaussLegendreQuadrature rule = new GaussLegendreQuadrature(n); + for (int k = 0; k <= 2*n - 1; k++) { + final int degree = k; + assertEquals(1.0 / (k + 1), rule.integrate((x) -> Math.pow(x, degree), 0, 1), 1E-14, + () -> "x^" + degree + " with " + rule.getNumPoints() + " points."); + } + } + final GaussLegendreQuadrature rule = new GaussLegendreQuadrature(2); + assertTrue(Math.abs(rule.integrate((x) -> x*x*x*x, 0, 1) - 0.2) > 1E-3, + "A rule of 2 points is not exact for a polynomial of degree 4."); + } + + /** + * A smooth function is integrated to the floating point accuracy with a few points. + */ + @Test + public void testSmoothFunction() { + final GaussLegendreQuadrature rule = new GaussLegendreQuadrature(10); + assertEquals(2, rule.integrate(Math::sin, 0, Math.PI), 1E-14); + assertEquals(Math.E - 1, rule.integrate(Math::exp, 0, 1), 1E-14); + } + + /** + * Exchanging the bounds changes the sign of the integral, and an empty interval gives zero. + */ + @Test + public void testBounds() { + final GaussLegendreQuadrature rule = new GaussLegendreQuadrature(10); + assertEquals(-2, rule.integrate(Math::sin, Math.PI, 0), 1E-14); + assertEquals( 0, rule.integrate(Math::sin, 1, 1), 0); + } + + /** + * The square root has an unbounded derivative at zero, so a single application of the rule is + * inaccurate. The adaptive integration halves the intervals near zero only, until the result is + * accurate. The sub-intervals reported to the consumer cover the interval without gap or overlap. + */ + @Test + public void testAdaptive() { + final GaussLegendreQuadrature rule = new GaussLegendreQuadrature(10); + final double expected = 2.0 / 3; + assertTrue(Math.abs(rule.integrate(Math::sqrt, 0, 1) - expected) > 1E-8, + "A single application of the rule is not accurate near the singularity."); + + final List<double[]> intervals = new ArrayList<>(); + final double result = rule.integrate(Math::sqrt, 0, 1, 1E-13, 40, + (start, end, integral) -> intervals.add(new double[] {start, end, integral})); + assertEquals(expected, result, 1E-12); + + double previousEnd = 0, sum = 0; + for (final double[] interval : intervals) { + assertEquals(previousEnd, interval[0], 0, "Sub-intervals shall be contiguous and in order."); + assertTrue(interval[1] > interval[0], "Sub-intervals shall not be empty."); + previousEnd = interval[1]; + sum += interval[2]; + } + assertEquals(1, previousEnd, 0, "Sub-intervals shall cover the whole interval."); + assertEquals(result, sum, 1E-15, "The integrals of the sub-intervals shall sum to the result."); + assertTrue(intervals.size() < 200, "Only the intervals near the singularity shall be halved."); + } + + /** + * Without halving, the adaptive integration integrates the two halves of the interval. + * The consumer is optional. + */ + @Test + public void testAdaptiveWithoutHalving() { + final GaussLegendreQuadrature rule = new GaussLegendreQuadrature(4); + final double halves = rule.integrate(Math::sqrt, 0, 0.5) + rule.integrate(Math::sqrt, 0.5, 1); + assertEquals(halves, rule.integrate(Math::sqrt, 0, 1, 1E-13, 0, null), 0); + } +}
