Merge pull request #27 from StaticBeagle/runge-kutta

Runge kutta

src/main/java/com/wildbitsfoundry/etk4j/control/TransferFunction.java

@@ -126,7 +126,7 @@ public Polynomial getNumerator() {
 
     /**
      * Numerator coefficients.
-     * @return An array containing the coefficients of the numerator in descending order.This is equivalent to
+     * @return An array containing the coefficients of the numerator in descending order.This is similar to
      * calling {@link #getNumerator()} and subsequently calling {@link Polynomial#getCoefficients()}.
      * However, this method pads the resulting array with zeros if the degree of the denominator is greater than the degree
      * of the numerator e.g.
@@ -135,7 +135,7 @@ public Polynomial getNumerator() {
      *     denominator: x^2 + 2 * x + 1 -> [1, 2, 1];
      * </pre>
      */
-    public double[] getNumeratorCoefficients() {
+    public double[] getNumeratorCoefficientsPadded() {
         Polynomial num = rf.getNumerator();
         int numOrder = num.degree();
         int denOrder = rf.getDenominator().degree();
@@ -148,6 +148,15 @@ public double[] getNumeratorCoefficients() {
         return result;
     }
 
+    /**
+     * Numerator coefficients.
+     * @return An array containing the coefficients of the numerator in descending order.This is equivalent to
+     * calling {@link #getNumerator()} and subsequently calling {@link Polynomial#getCoefficients()}.
+     */
+    public double[] getNumeratorCoefficients() {
+        return this.rf.getNumerator().getCoefficients();
+    }
+
     /**
      * Denominator {@link Polynomial}.
      * @return The denominator polynomial.
@@ -160,14 +169,23 @@ public Polynomial getDenominator() {
      * Denominator coefficients.
      * @return An array containing the coefficients of the denominator in descending order. This is equivalent to
      * calling {@link #getDenominator()} and subsequently calling {@link Polynomial#getCoefficients()}.
+     */
+    public double[] getDenominatorCoefficients() {
+        return  rf.getDenominator().getCoefficients();
+    }
+
+    /**
+     * Denominator coefficients.
+     * @return An array containing the coefficients of the denominator in descending order. This is similar to
+     * calling {@link #getDenominator()} and subsequently calling {@link Polynomial#getCoefficients()}.
      * However, this method pads the resulting array with zeros if the degree of the numerator is greater than the degree
      * of the denominator e.g.
      * <pre>
      *     numerator: x^2 + 2 * x + 1 -&gt; [1, 2, 1];<br>
      *     denominator: x + 1         -&gt; [0, 1, 1];
      * </pre>
      */
-    public double[] getDenominatorCoefficients() {
+    public double[] getDenominatorCoefficientsPadded() {
         Polynomial den = rf.getDenominator();
         int numOrder = rf.getNumerator().degree();
         int denOrder = den.degree();

src/main/java/com/wildbitsfoundry/etk4j/math/MathETK.java

@@ -174,7 +174,7 @@ public static class FRexpResult {
 	 * @return The mantissa and exponent of a number such that number = m * 2^e.
 	 * @see <a href="https://stackoverflow.com/a/3946294/6383857">https://stackoverflow.com/a/3946294/6383857</a>
 	 */
-	public static FRexpResult frexp(double value) {
+	public static FRexpResult frExp(double value) {
 		FRexpResult result = new FRexpResult();
 
 		result.exponent = 0;
@@ -229,4 +229,26 @@ public static FRexpResult frexp(double value) {
 
 		return result;
 	}
+
+	/**
+	 * Combinations of picking k unordered outcomes from n possibilities
+	 * @param n the number of things
+	 * @param k the number of elements taken
+	 * @return The total number of combinations
+	 */
+	public static long combinations(long n, long k) {
+		return factorial(n) / (factorial(k) * factorial(n - k));
+	}
+
+	/**
+	 * Factorial of a number
+	 * @param n the number
+	 * @return {@code n!}
+	 */
+	public static long factorial(long n) {
+		if(n == 0) {
+			return 1;
+		}
+		return n * factorial(n - 1);
+	}
 }

src/main/java/com/wildbitsfoundry/etk4j/math/calculus/odesolvers/ODESystemOfEquations.java

@@ -0,0 +1,5 @@
+package com.wildbitsfoundry.etk4j.math.calculus.odesolvers;
+
+public interface ODESystemOfEquations {
+    double[] evaluateAt(double t, double[] y);
+}

src/main/java/com/wildbitsfoundry/etk4j/math/calculus/odesolvers/OdeSolver.java

@@ -0,0 +1,86 @@
+package com.wildbitsfoundry.etk4j.math.calculus.odesolvers;
+
+import com.wildbitsfoundry.etk4j.util.Tuples;
+
+import java.util.Arrays;
+// TODO document this class
+public abstract class OdeSolver {
+
+    protected double t;
+    protected Double tOld;
+    protected double[] y;
+    protected ODESystemOfEquations systemOfEquations;
+    protected double tBound;
+    protected double direction;
+    protected String status; // TODO change to enum
+    protected int n;
+
+    public OdeSolver(ODESystemOfEquations systemOfEquations, double t0, double[] y0, Double tBound) {
+        this.t = t0;
+        if (!Arrays.stream(y0).allMatch(Double::isFinite)) {
+            throw new IllegalArgumentException("All components of the initial state y0 must be finite.");
+        }
+        this.systemOfEquations = systemOfEquations;
+        this.direction = (int) Math.signum(tBound - t0);
+        this.y = y0;
+        this.n = y0.length;
+        this.status = "running";
+        // nfev
+        // njev
+        // nlu
+    }
+
+    private Double stepSize() {
+        return tOld == null ? null : Math.abs(t - tOld);
+    }
+
+    public String step() {
+        if (this.status != "running") {
+            throw new RuntimeException("Attempt to step on a failed or finished solver.");
+        }
+
+        if (this.n == 0 || this.t == this.tBound) {
+            this.tOld = this.t;
+            this.t = this.tBound;
+            this.status = "finished";
+            return null;
+        }
+        Tuples.Tuple2<Boolean, String> result = stepImpl();
+
+        if (!result.getItem1()) {
+            this.status = "failed";
+        } else {
+            this.tOld = this.t;
+            if (this.direction * (this.t - this.tBound) >= 0) {
+                status = "finished";
+            }
+        }
+        return result.getItem2();
+    }
+
+    protected abstract Tuples.Tuple2<Boolean, String> stepImpl();
+
+    // TODO dense output
+//    public Object getDenseOutput() {
+//        if (tOld == null) {
+//            throw new RuntimeException("Dense output is available after a successful step was made.");
+//        }
+//
+//        if (n == 0 || t == tOld) {
+//            // return ConstantDenseOutput(told, t, y);
+//            return null; // TODO
+//        } else {
+//            return getDenseOutputImpl();
+//        }
+//    }
+//
+//    protected abstract DenseOutput getDenseOutputImpl();
+//
+//    public static class DenseOutput {
+//
+//    }
+//
+//    public static class ConstantDenseOutput extends DenseOutput {
+//
+//    }
+}