Merge pull request #25 from StaticBeagle/romberg-integral

staticbeagle · web-flow · commit 968f9de7973b · 2024-12-03T21:23:18.000-06:00
Added Romberg integrals, Jacobi solver, added release profile and clean-up
diff --git a/.github/workflows/release.yml b/.github/workflows/release.yml
@@ -26,4 +26,4 @@ jobs:
           gpg-private-key: ${{ secrets.GPG_PRIVATE_KEY }}
           gpg-passphrase: GPG_PASSPHRASE
       - name: Publish package
-        run: mvn -B deploy
+        run: mvn -B deploy -P release
diff --git a/README.md b/README.md
@@ -2,7 +2,6 @@
 
 [![Java CI with Maven](https://github.com/StaticBeagle/ETK4J/actions/workflows/maven.yml/badge.svg)](https://github.com/StaticBeagle/ETK4J/actions/workflows/maven.yml)
 
----------------------------------------------------------------------------
 ## Engineering Toolkit for Java
 
 The purpose of this project is to create a library that can be used to prototype solutions to engineering problems. One 
@@ -27,15 +26,22 @@ can be found in:
 
     src/main/java
 
----------------------------------------------------------------------------
 ## Maven Central
 
-ETK4J can be included from Maven Central
+ETK4J can be included from Maven Central.
+
+Maven
 
 ```xml
-<groupId>com.wildbitsfoundry</groupId>
-<artifactId>etk4j</artifactId>
-<version>2.0.0</version>
+<dependency>
+    <groupId>com.wildbitsfoundry</groupId>
+    <artifactId>etk4j</artifactId>
+    <version>2.0.0</version>
+</dependency>
+```
+Gradle
+```bash
+implementation 'com.wildbitsfoundry:etk4j:2.0.0'
 ```
 
 Requirements
diff --git a/pom.xml b/pom.xml
@@ -2,7 +2,7 @@
     <modelVersion>4.0.0</modelVersion>
     <groupId>com.wildbitsfoundry</groupId>
     <artifactId>etk4j</artifactId>
-    <version>2.0.0</version>
+    <version>2.1.0</version>
     <name>Engineering Toolkit for Java</name>
     <description>Tools and implementation of mathematical methods or backing math for engineering problems and
         applications.
@@ -12,7 +12,7 @@
         <project.build.sourceEncoding>UTF-8</project.build.sourceEncoding>
     </properties>
     <url>
-        https://github.com/StaticBeagle/ETK4J
+        https://wildbitsfoundry.com/etk4j
     </url>
 
     <developers>
@@ -32,9 +32,8 @@
     </licenses>
 
     <scm>
-        <url>
-            scm:git:git@github.com:StaticBeagle/ETK4J.git
-        </url>
+        <connection>scm:git:https://github.com/StaticBeagle/ETK4J.git</connection>
+        <url>https://github.com/StaticBeagle/ETK4J</url>
     </scm>
 
     <dependencies>
@@ -47,6 +46,30 @@
         </dependency>
     </dependencies>
 
+    <profiles>
+        <profile>
+            <id>release</id>
+            <build>
+                <plugins>
+                    <plugin>
+                        <groupId>org.apache.maven.plugins</groupId>
+                        <artifactId>maven-gpg-plugin</artifactId>
+                        <version>3.2.7</version>
+                        <executions>
+                            <execution>
+                                <id>sign-artifacts</id>
+                                <phase>verify</phase>
+                                <goals>
+                                    <goal>sign</goal>
+                                </goals>
+                            </execution>
+                        </executions>
+                    </plugin>
+                </plugins>
+            </build>
+        </profile>
+    </profiles>
+
     <build>
         <plugins>
             <plugin>
@@ -100,20 +123,6 @@
                     </execution>
                 </executions>
             </plugin>
-            <plugin>
-                <groupId>org.apache.maven.plugins</groupId>
-                <artifactId>maven-gpg-plugin</artifactId>
-                <version>3.2.7</version>
-                <executions>
-                    <execution>
-                        <id>sign-artifacts</id>
-                        <phase>verify</phase>
-                        <goals>
-                            <goal>sign</goal>
-                        </goals>
-                    </execution>
-                </executions>
-            </plugin>
         </plugins>
     </build>
 </project>
diff --git a/src/main/java/com/wildbitsfoundry/etk4j/math/calculus/Integrals.java b/src/main/java/com/wildbitsfoundry/etk4j/math/calculus/Integrals.java
@@ -119,6 +119,56 @@ public static double simpson(UnivariateFunction func, double a, double b, int n)
         return simpson((x, o) -> func.evaluateAt(x), a, b, n);
     }
 
+    /**
+     * Computes the approximate definite integral from a to b using the Romberg and Richardson extrapolation.
+     *
+     * @param func   The function to be integrated.
+     * @param a      The lower limit of the integral.
+     * @param b      The upper limit of the integral.
+     * @param maxIterations The maximum number of iterations
+     * @param params Optional parameters passed to {@code func}.
+     * @return The approximate definite integral of {@code func} from a to b.
+     */
+    public static double romberg(BiFunction<Double, Object[], Double> func, double a, double b, int maxIterations, Object... params) {
+        double[][] R = new double[maxIterations + 1][maxIterations + 1];
+
+        R[0][0] = trapz(func, a, b, 1, params);
+        for(int i = 0; i <= maxIterations; i++) {
+            int n = (int) Math.pow(2, i);
+            R[i][0] = trapz(func, a, b, n, params);
+
+            for(int j = 1; j <= i; j++) {
+                R[i][j] = (Math.pow(4, j) * R[i][j - 1] - R[i - 1][j - 1]) / (Math.pow(4, j) - 1);
+            }
+        }
+        return R[maxIterations][maxIterations];
+    }
+
+    /**
+     * Computes the approximate definite integral from a to b using the Romberg and Richardson extrapolation.
+     *
+     * @param func   The function to be integrated.
+     * @param a      The lower limit of the integral.
+     * @param b      The upper limit of the integral.
+     * @param maxIterations The maximum number of iterations
+     * @return The approximate definite integral of {@code func} from a to b.
+     */
+    public static double romberg(UnivariateFunction func, double a, double b, int maxIterations) {
+        return romberg((x, o) -> func.evaluateAt(x), a, b, maxIterations);
+    }
+
+    /**
+     * Computes the approximate definite integral from a to b using the Romberg and Richardson extrapolation.
+     *
+     * @param func   The function to be integrated.
+     * @param a      The lower limit of the integral.
+     * @param b      The upper limit of the integral.
+     * @return The approximate definite integral of {@code func} from a to b.
+     */
+    public static double romberg(UnivariateFunction func, double a, double b) {
+        return romberg((x, o) -> func.evaluateAt(x), a, b, 10);
+    }
+
     /**
      * Computes the approximate definite integral from a to b using 5 point Gaussian quadrature.
      *
diff --git a/src/main/java/com/wildbitsfoundry/etk4j/math/interpolation/QuadraticSpline.java b/src/main/java/com/wildbitsfoundry/etk4j/math/interpolation/QuadraticSpline.java
@@ -1,6 +1,6 @@
 package com.wildbitsfoundry.etk4j.math.interpolation;
 
-import com.wildbitsfoundry.etk4j.math.linearalgebra.GaussianElimination;
+import com.wildbitsfoundry.etk4j.math.linearalgebra.GaussianEliminationSolver;
 import com.wildbitsfoundry.etk4j.math.linearalgebra.MatrixDense;
 
 import java.util.Arrays;
@@ -41,7 +41,7 @@ public static QuadraticSpline newNaturalSplineInPlace(double[] x, double[] y) {
         A.set(2 * n - 4, n - 1, 1); // c[0] = 0
         A.set(2 * n - 3, 2 * n - 3, 1); // c[n - 1] = 0
 
-        double[] solution = GaussianElimination.solve(A, b);
+        double[] solution = GaussianEliminationSolver.solve(A, b);
         // compute coefficients
         double[] coefficients = new double[(n - 1) * 3];
         for (int i = 0, j = 0; i < n - 1; ++i, ++j) {
@@ -83,7 +83,7 @@ public static QuadraticSpline newClampedSplineInPlace(double[] x, double[] y, do
         A.set(2 * n - 3, 2 * n - 3, 2 * (x[n - 1] - x[n - 2]));
         b[2 * n - 3] = mn; // derivative at the end
 
-        double[] solution = GaussianElimination.solve(A, b);
+        double[] solution = GaussianEliminationSolver.solve(A, b);
         // compute coefficients
         double[] coefficients = new double[(n - 1) * 3];
         for (int i = 0, j = 0; i < n - 1; ++i, ++j) {
diff --git a/src/main/java/com/wildbitsfoundry/etk4j/math/linearalgebra/GaussSeidel.java b/src/main/java/com/wildbitsfoundry/etk4j/math/linearalgebra/GaussSeidel.java
diff --git a/src/main/java/com/wildbitsfoundry/etk4j/math/linearalgebra/GaussSeidelSolver.java b/src/main/java/com/wildbitsfoundry/etk4j/math/linearalgebra/GaussSeidelSolver.java
@@ -0,0 +1,8 @@
+package com.wildbitsfoundry.etk4j.math.linearalgebra;
+// TODO javadocs
+public class GaussSeidelSolver extends SuccessiveOverRelaxationSolver {
+
+    public GaussSeidelSolver(MatrixDense A, double[] b) {
+        super(A, b, 1.0);
+    }
+}
diff --git a/src/main/java/com/wildbitsfoundry/etk4j/math/linearalgebra/GaussianEliminationSolver.java b/src/main/java/com/wildbitsfoundry/etk4j/math/linearalgebra/GaussianEliminationSolver.java
@@ -1,8 +1,8 @@
 package com.wildbitsfoundry.etk4j.math.linearalgebra;
+// TODO javadocs
+public final class GaussianEliminationSolver {
 
-public final class GaussianElimination {
-
-    private GaussianElimination() {}
+    private GaussianEliminationSolver() {}
 
     public static double[] solve(MatrixDense A, double[] b) {
         return solve(A, b, 1e-10);
diff --git a/src/main/java/com/wildbitsfoundry/etk4j/math/linearalgebra/JacobiMethodSolver.java b/src/main/java/com/wildbitsfoundry/etk4j/math/linearalgebra/JacobiMethodSolver.java
@@ -0,0 +1,84 @@
+package com.wildbitsfoundry.etk4j.math.linearalgebra;
+
+import com.wildbitsfoundry.etk4j.util.DoubleArrays;
+// TODO javadocs
+public class JacobiMethodSolver {
+
+    private double[] b;
+    private MatrixDense A;
+    private int iterationLimit = 100;
+    private double tol = 1e-9;
+    private double[] x0 = null;
+
+    public JacobiMethodSolver(MatrixDense A, double[] b) {
+        this.A = A;
+        this.b = b;
+    }
+
+    public JacobiMethodSolver iterationLimit(int iterationLimit) {
+        this.iterationLimit = iterationLimit;
+        return this;
+    }
+
+    public JacobiMethodSolver tolerance(double tolerance) {
+        this.tol = tolerance;
+        return this;
+    }
+
+    public JacobiMethodSolver initialGuess(double[] x0) {
+        this.x0 = x0;
+        return this;
+    }
+
+    public IterativeSolverResults<double[]> solve() {
+        double[] x = x0 == null ? new double[b.length] : x0;
+        int n = A.getRowCount();
+        double[] xNew = new double[n];
+        int k;
+        double error = Double.NaN;
+        for (k = 0; k < iterationLimit; k++) {
+            for (int i = 0; i < n; i++) {
+                double sum = 0;
+                for (int j = 0; j < n; j++) {
+                    if (j != i) {
+                        sum += A.unsafeGet(i, j) * x[j];
+                    }
+                }
+                xNew[i] = (b[i] - sum) / A.unsafeGet(i, i);
+            }
+            // Check for convergence
+            error = DoubleArrays.norm2(DoubleArrays.subtractElementWise(xNew, x));
+            if(error < tol) {
+                IterativeSolverResults<double[]> solverResults = new IterativeSolverResults<>();
+                solverResults.setSolverStatus("Converged");
+                solverResults.setHasConverged(true);
+                solverResults.setError(error);
+                solverResults.setValue(xNew);
+                solverResults.setNumberOfIterations(k + 1);
+                return solverResults;
+            }
+            System.arraycopy(xNew, 0, x, 0, xNew.length);
+        }
+        IterativeSolverResults<double[]> solverResults = new IterativeSolverResults<>();
+        solverResults.setSolverStatus("Maximum number of iterations exceeded");
+        solverResults.setHasConverged(false);
+        solverResults.setError(error);
+        solverResults.setValue(x);
+        solverResults.setNumberOfIterations(k + 1);
+        return solverResults;
+
+    }
+    public static void main(String[] args) {
+        double[][] A = {{2, -1, 0}, {-1, 2, -1}, {0, -1, 2}};
+        double[] b = {1, 0, 1};
+        double[] x = {0, 0, 0}; // Initial guess
+        int maxIterations = 100;
+        double tolerance = 0.0001;
+        double[] solution = new JacobiMethodSolver(new MatrixDense(A), b)
+                .initialGuess(x).iterationLimit(maxIterations).tolerance(tolerance).solve().getValue();
+        System.out.println("Solution:");
+        for (double val : solution) {
+            System.out.println(val);
+        }
+    }
+}
diff --git a/src/main/java/com/wildbitsfoundry/etk4j/math/linearalgebra/SuccessiveOverRelaxationSolver.java b/src/main/java/com/wildbitsfoundry/etk4j/math/linearalgebra/SuccessiveOverRelaxationSolver.java
@@ -1,32 +1,31 @@
 package com.wildbitsfoundry.etk4j.math.linearalgebra;
 
-import java.util.Arrays;
 // TODO javadocs
-public class SuccessiveOverRelaxation {
+public class SuccessiveOverRelaxationSolver {
     private double w;
     private double[] b;
     private MatrixDense A;
     private int iterationLimit = 100;
     private double tol = 1e-9;
     private double[] x0 = null;
 
-    public SuccessiveOverRelaxation(MatrixDense A, double[] b, double w) {
+    public SuccessiveOverRelaxationSolver(MatrixDense A, double[] b, double w) {
         this.A = A;
         this.b = b;
         this.w = w;
     }
 
-    public SuccessiveOverRelaxation iterationLimit(int iterationLimit) {
+    public SuccessiveOverRelaxationSolver iterationLimit(int iterationLimit) {
         this.iterationLimit = iterationLimit;
         return this;
     }
 
-    public SuccessiveOverRelaxation tolerance(double tolerance) {
+    public SuccessiveOverRelaxationSolver tolerance(double tolerance) {
         this.tol = tolerance;
         return this;
     }
 
-    public SuccessiveOverRelaxation initialConditions(double[] x0) {
+    public SuccessiveOverRelaxationSolver initialGuess(double[] x0) {
         this.x0 = x0;
         return this;
     }
diff --git a/src/test/java/com/wildbitsfoundry/etk4j/math/calculus/IntegralsTest.java b/src/test/java/com/wildbitsfoundry/etk4j/math/calculus/IntegralsTest.java
@@ -35,6 +35,13 @@ public void testSimpson() {
 		assertEquals(0.8281163288433171, simp, 1e-12);
 	}
 
+	@Test
+	public void testRomberg() {
+		UnivariateFunction fx = x -> Math.sin(x * x);
+		double romb = Integrals.romberg(fx, 0, Math.PI / 2.0);
+		assertEquals(0.8281163288428953, romb, 1e-12);
+	}
+
 	@Test
 	public void testGaussianQuadrature() {
 		UnivariateFunction fx = x -> Math.sin(x * x);
@@ -45,7 +52,7 @@ public void testGaussianQuadrature() {
 		qadrat = Integrals.gaussianQuadrature(fx, Math.PI / 4.0, Math.PI / 2.0);
 		assertEquals(0.5832680904501192, qadrat, 1e-12);
 	}
-	
+
 	@Test
 	public void testAdaptiveGaussianQuadrature() {
 		double[] e = new double[2];
diff --git a/src/test/java/com/wildbitsfoundry/etk4j/math/linearalgebra/GaussSeidelSolverTest.java b/src/test/java/com/wildbitsfoundry/etk4j/math/linearalgebra/GaussSeidelSolverTest.java
diff --git a/src/test/java/com/wildbitsfoundry/etk4j/math/linearalgebra/JacobiMethodSolverTest.java b/src/test/java/com/wildbitsfoundry/etk4j/math/linearalgebra/JacobiMethodSolverTest.java
