diff --git a/.classpath b/.classpath deleted file mode 100644 index bdb2412f..00000000 --- a/.classpath +++ /dev/null @@ -1,6 +0,0 @@ - - - - - - diff --git a/src/eva2/tools/math/Mathematics.java b/src/eva2/tools/math/Mathematics.java index d0f7549e..6c6bd3f5 100644 --- a/src/eva2/tools/math/Mathematics.java +++ b/src/eva2/tools/math/Mathematics.java @@ -18,23 +18,25 @@ import eva2.tools.math.interpolation.SplineInterpolation; public class Mathematics { /** * Computes the full adjoint matrix. - * + * * @param a * @return */ public static double[][] adjoint(double[][] a) { - if (a == null) return null; - if (a.length != a[0].length) return null; + if (a == null) + return null; + if (a.length != a[0].length) + return null; double[][] b = new double[a.length][a.length]; for (int i = 0; i < a.length; i++) for (int j = 0; j < a.length; j++) b[i][j] = adjoint(a, i, j); return b; } - + /** * Computes the adjoint of the matrix element at the position (k, l). - * + * * @param a * @param k * @param l @@ -43,53 +45,58 @@ public class Mathematics { public static double adjoint(double[][] a, int k, int l) { return Math.pow(-1, k + l + 2) * determinant(submatrix(a, k, l)); } - + /** * This computes the determinant of the given matrix - * + * * @param matrix * @return The determinant or null if there is no determinant (if the matrix * is not square). */ public static double determinant(double[][] matrix) { - if (matrix == null) return 0; - if (matrix.length != matrix[0].length) return 0; - if (matrix.length == 1) return matrix[0][0]; + if (matrix == null) + return 0; + if (matrix.length != matrix[0].length) + return 0; + if (matrix.length == 1) + return matrix[0][0]; if (matrix.length == 2) return matrix[0][0] * matrix[1][1] - matrix[0][1] * matrix[1][0]; if (matrix.length == 3) return matrix[0][0] * matrix[1][1] * matrix[2][2] + matrix[0][1] - * matrix[1][2] * matrix[2][0] + matrix[0][2] * matrix[1][0] - * matrix[2][1] - matrix[2][0] * matrix[1][1] * matrix[0][2] - - matrix[2][1] * matrix[1][2] * matrix[0][0] - matrix[2][2] - * matrix[1][0] * matrix[0][1]; + * matrix[1][2] * matrix[2][0] + matrix[0][2] * matrix[1][0] + * matrix[2][1] - matrix[2][0] * matrix[1][1] * matrix[0][2] + - matrix[2][1] * matrix[1][2] * matrix[0][0] - matrix[2][2] + * matrix[1][0] * matrix[0][1]; double det = 0; for (int k = 0; k < matrix.length; k++) { - if (matrix[0][k] != 0) det += matrix[0][k] * adjoint(matrix, 0, k); + if (matrix[0][k] != 0) + det += matrix[0][k] * adjoint(matrix, 0, k); } return det; } - + /** * Computes the root-Distance function. For example root = 2 gives the * Euclidian Distance. - * + * * @param x - * a vector + * a vector * @param y - * another vector + * another vector * @param root - * what kind of distance funktion + * what kind of distance funktion * @return the distance of x and y * @throws Exception - * if x and y have different dimensions an exception is - * thrown. + * if x and y have different dimensions an exception is thrown. */ public static double dist(double[] x, double[] y, int root) { if (x.length != y.length) - throw new RuntimeException("The vectors x and y must have the same dimension"); - if (root == 0) throw new RuntimeException("There is no 0-root!"); + throw new RuntimeException( + "The vectors x and y must have the same dimension"); + if (root == 0) + throw new RuntimeException("There is no 0-root!"); double d = 0; for (int i = 0; i < x.length; i++) d += Math.pow(Math.abs(x[i] - y[i]), root); @@ -98,30 +105,30 @@ public class Mathematics { /** * Computes the euclidian distance function. - * + * * @param x - * a vector + * a vector * @param y - * another vector + * another vector * @param root - * what kind of distance funktion + * what kind of distance funktion * @return the distance of x and y * @throws Exception - * if x and y have different dimensions an exception is - * thrown. + * if x and y have different dimensions an exception is thrown. */ - public static double euclidianDist(double[] x, double[] y) { + public static double euclidianDist(double[] x, double[] y) { if (x.length != y.length) - throw new RuntimeException("The vectors x and y must have the same dimension"); + throw new RuntimeException( + "The vectors x and y must have the same dimension"); double d = 0; for (int i = 0; i < x.length; i++) d += Math.pow(Math.abs(x[i] - y[i]), 2); return Math.sqrt(d); } -/** - * Expand a vector to a higher dimension (len) by filling it up - * with a constant value. + /** + * Expand a vector to a higher dimension (len) by filling it up with a + * constant value. * * @param x * @param len @@ -129,15 +136,17 @@ public class Mathematics { * @return */ public static double[] expandVector(double[] x, int len, double v) { - if (len <= x.length ) {// not really an error, just perform identity -// System.err.println("Error, invalid length in expandVector, expecting l>" + x.length); + if (len <= x.length) {// not really an error, just perform identity + // System.err.println("Error, invalid length in expandVector, expecting l>" + // + x.length); return x; - } else { - double[] expanded = new double[len]; - System.arraycopy(x, 0, expanded, 0, x.length); - for (int i=x.length; iupper)) return false; + if (v < lower || (v > upper)) + return false; return true; } @@ -358,32 +389,37 @@ public class Mathematics { * * @param x * @param range - * @return true if the vector lies within the range, else false + * @return true if the vector lies within the range, else false */ public static boolean isInRange(double[] x, double[][] range) { - for (int i=0; irange[i][1])) return false; + for (int i = 0; i < x.length; i++) { + if (x[i] < range[i][0] || (x[i] > range[i][1])) + return false; } return true; } /** - * Returns false if a vector contains NaN, its squared sum is NaN - * or the absolute sum is smaller than 10^-18. + * Returns false if a vector contains NaN, its squared sum is NaN or the + * absolute sum is smaller than 10^-18. + * * @param d * @return */ public static boolean isValidVec(double[] d) { double sum = 0; for (int i = 0; i < d.length; i++) { - if (Double.isNaN(d[i])) return false; - sum += Math.pow(d[i],2); + if (Double.isNaN(d[i])) + return false; + sum += Math.pow(d[i], 2); } - if (Double.isNaN(sum)) return false; - if (Math.abs(sum) < 0.000000000000000001) return false; + if (Double.isNaN(sum)) + return false; + if (Math.abs(sum) < 0.000000000000000001) + return false; return true; } - + /** * @param f0 * @param f1 @@ -397,46 +433,52 @@ public class Mathematics { /** * This method gives a linear interpolation of the function values of the * given argument/function value pairs. - * + * * @param x - * The argument at the point with unknown function value + * The argument at the point with unknown function value * @param x0 - * The argument at the last position with a function value + * The argument at the last position with a function value * @param x1 - * The argument at the next known fuction value + * The argument at the next known fuction value * @param f0 - * The function value at the position x0 + * The function value at the position x0 * @param f1 - * The function value at the position x1 + * The function value at the position x1 * @return The function value at position x given by linear interpolation. */ public static double linearInterpolation(double x, double x0, double x1, double f0, double f1) { - if (x1 == x0) return f0; + if (x1 == x0) + return f0; return lerp(f0, f1, (x - x0) / (x1 - x0)); } - + public static double max(double[] vals) { double maxVal = vals[0]; - for (int i=1; i dblArrList, boolean interpolate) { - java.util.Collections.sort(dblArrList, new DoubleArrayComparator()); // by default, the comparator uses pareto dominance - + java.util.Collections.sort(dblArrList, new DoubleArrayComparator()); // by + // default, + // the + // comparator + // uses + // pareto + // dominance + int len = dblArrList.size(); - if (len % 2 != 0) return dblArrList.get((len-1) / 2); + if (len % 2 != 0) + return dblArrList.get((len - 1) / 2); else { - double[] med = dblArrList.get(len/2).clone(); + double[] med = dblArrList.get(len / 2).clone(); if (interpolate) { - vvAdd(med, dblArrList.get((len/2)+1), med); + vvAdd(med, dblArrList.get((len / 2) + 1), med); svDiv(2, med, med); } return med; @@ -501,45 +560,52 @@ public class Mathematics { public static double min(double[] vals) { double minVal = vals[0]; - for (int i=1; irange.length) System.err.println("Invalid vector length, x is longer than range! (Mathematics.projectToRange)"); - for (int i=0; i range.length) + System.err + .println("Invalid vector length, x is longer than range! (Mathematics.projectToRange)"); + for (int i = 0; i < x.length; i++) { + if (x[i] < range[i][0]) { viols++; - x[i]=range[i][0]; - } else if (x[i]>range[i][1]) { + x[i] = range[i][0]; + } else if (x[i] > range[i][1]) { viols++; - x[i]=range[i][1]; + x[i] = range[i][1]; } } return viols; @@ -628,64 +699,72 @@ public class Mathematics { * @return the closest value to v within [min,max] */ public static double projectValue(double v, double min, double max) { - if (vmax) { + } else if (v > max) { return max; - } else return v; + } else + return v; } - + /** - * Create a random vector, the components will be set to gaussian distributed - * values with mean zero and the given standard deviation. - * - * @param dim the desired dimension - * @param stdDev the gaussian standard deviation - * @return random vector + * Create a random vector, the components will be set to gaussian + * distributed values with mean zero and the given standard deviation. + * + * @param dim + * the desired dimension + * @param stdDev + * the gaussian standard deviation + * @return random vector */ public static double[] randomVector(int dim, double stdDev) { - double[] vect = new double[dim]; + double[] vect = new double[dim]; for (int j = 0; j < vect.length; j++) { - vect[j] = RNG.gaussianDouble(stdDev); + vect[j] = RNG.gaussianDouble(stdDev); } return vect; } - + /** - * Reflect the entries of x which violate the bounds to within the range. - * Return the number of violating dimensions. - * - * @param x - * @param range - * @return the number of violating dimensions - */ + * Reflect the entries of x which violate the bounds to within the range. + * Return the number of violating dimensions. + * + * @param x + * @param range + * @return the number of violating dimensions + */ public static int reflectBounds(double[] x, double[][] range) { - int viols=0; + int viols = 0; double d = 0.; - for (int i=0; i dimLen) d -= dimLen; // avoid violating the other bound immediately - x[i]=range[i][0]+d; - } else if (x[i]>range[i][1]) { + d = range[i][0] - x[i]; + while (d > dimLen) + d -= dimLen; // avoid violating the other bound + // immediately + x[i] = range[i][0] + d; + } else if (x[i] > range[i][1]) { viols++; - d = x[i]-range[i][1]; - while (d>dimLen) d -= dimLen; // avoid violating the other bound immediately - x[i]=range[i][1]-d; + d = x[i] - range[i][1]; + while (d > dimLen) + d -= dimLen; // avoid violating the other bound + // immediately + x[i] = range[i][1] - d; } } } return viols; } - + /** - * Simple version of reflection of a value moving by a step and bouncing - * of min and max values like a pool ball. Precondition is min <= val <= max, + * Simple version of reflection of a value moving by a step and bouncing of + * min and max values like a pool ball. Precondition is min <= val <= max, * post condition is min <= retVal <= max. * * @param val @@ -694,61 +773,66 @@ public class Mathematics { * @param max * @return */ - public static double reflectValue(double val, double step, double min, double max) { - while (step > (max-min)) step -= (max-min); - if ((val + step) > max) + public static double reflectValue(double val, double step, double min, + double max) { + while (step > (max - min)) + step -= (max - min); + if ((val + step) > max) return (2 * max - val - step); - if ((val + step) < min) - return (2 * min - val - step); + if ((val + step) < min) + return (2 * min - val - step); return (val += step); } - + /** * Computes the relative distance of vector x to vector y. Therefore the * difference of x[i] and y[i] is divided by y[i] for every i. If y[i] is * zero, the default value def is used instead. The sum of these differences * gives the distance function. - * + * * @param x - * A vector + * A vector * @param y - * The reference vector + * The reference vector * @param def - * The default value to be use to avoid division by zero. + * The default value to be use to avoid division by zero. * @return The relative distance of x to y. * @throws Exception */ public static double relDist(double[] x, double[] y, double def) - throws Exception { + throws Exception { if (x.length != y.length) - throw new Exception("The vectors x and y must have the same dimension"); + throw new Exception( + "The vectors x and y must have the same dimension"); double d = 0; for (int i = 0; i < x.length; i++) if (y[i] != 0) d += Math.pow(((x[i] - y[i]) / y[i]), 2); - else d += def; + else + d += def; return d; } - + /** - * Reset single axis rotation matrix to unity. - */ - public static void resetRotationEntriesSingleAxis(Matrix tmp, int i, int j) { - tmp.set(i, i, 1); - tmp.set(i, j, 0); - tmp.set(j, i, 0); - tmp.set(j, j, 1); - } - - public static void revertArray(Object[] src, Object[] dst) { - if (dst.length>=src.length) { - for (int i=0; i= src.length) { + for (int i = 0; i < src.length; i++) { + dst[src.length - i - 1] = src[i]; + } + } else + System.err.println("Mismatching array lengths!"); + } + + /** * Rotate the vector by angle alpha around axis i/j * * @param vect @@ -759,77 +843,85 @@ public class Mathematics { public static void rotate(double[] vect, double alpha, int i, int j) { double xi = vect[i]; double xj = vect[j]; - vect[i] = (xi*Math.cos(alpha))-(xj*Math.sin(alpha)); - vect[j] = (xi*Math.sin(alpha))+(xj*Math.cos(alpha)); + vect[i] = (xi * Math.cos(alpha)) - (xj * Math.sin(alpha)); + vect[j] = (xi * Math.sin(alpha)) + (xj * Math.cos(alpha)); } - - /** - * Rotate a given double vector using a rotation matrix. If the matrix - * is null, x will be returned unchanged. Matrix dimensions must fit. - * - * @param x - * @param rotMatrix - * @return the rotated vector - */ - public static double[] rotate(double[] x, Matrix rotMatrix) { - if (rotMatrix!=null) { - Matrix resVec = rotMatrix.times(new Matrix(x, x.length)); - x = resVec.getColumnPackedCopy(); - return x; - } else return x; - } - /** + /** + * Rotate a given double vector using a rotation matrix. If the matrix is + * null, x will be returned unchanged. Matrix dimensions must fit. + * + * @param x + * @param rotMatrix + * @return the rotated vector + */ + public static double[] rotate(double[] x, Matrix rotMatrix) { + if (rotMatrix != null) { + Matrix resVec = rotMatrix.times(new Matrix(x, x.length)); + x = resVec.getColumnPackedCopy(); + return x; + } else + return x; + } + + /** * Rotate the vector along all axes by angle alpha or a uniform random value * in [-alpha, alpha] if randomize is true. - * + * * @param vect * @param alpha * @param randomize */ - public static void rotateAllAxes(double[] vect, double alpha, boolean randomize) { - for (int i=0; i1) or reduced (fact<1) by the defined ratio around the center. - * - * @param rangeScaleFact - * @param range - */ + * Scale a range by the given factor, meaning that the interval in each + * dimension is extended (fact>1) or reduced (fact<1) by the defined ratio + * around the center. + * + * @param rangeScaleFact + * @param range + */ public static void scaleRange(double rangeScaleFact, double[][] range) { - double[] intervalLengths=Mathematics.getAbsRange(range); - double[] tmpInts=Mathematics.svMult(rangeScaleFact, intervalLengths); - Mathematics.vvSub(tmpInts, intervalLengths, tmpInts); // this is what must be added to range interval - for (int i=0; i - DESystem */ public interface DESSolver extends Serializable { diff --git a/src/eva2/tools/math/des/RKSolver.java b/src/eva2/tools/math/des/RKSolver.java index 364544bc..fa8336a2 100644 --- a/src/eva2/tools/math/des/RKSolver.java +++ b/src/eva2/tools/math/des/RKSolver.java @@ -12,6 +12,7 @@ import eva2.tools.math.Mathematics; * @author Andreas Dräger * @author Marcel Kronfeld * @version 1.0 Status: works, but numerical inaccurate + * @depend - - Mathematics */ public class RKSolver implements DESSolver, Serializable { /** @@ -89,109 +90,6 @@ public class RKSolver implements DESSolver, Serializable { return unstableFlag; } - /** - * @param DES - * @param h - * @param x - * @param Ytemp - * @return - */ - public double[] rkTerm(DESystem DES, double h, double x, double[] Ytemp) { - double[][] K = new double[4][]; - K[0] = Mathematics.svMult(h, DES.getValue(x, Ytemp)); - K[1] = Mathematics.svMult(h, DES.getValue(x + h / 2, Mathematics.vvAdd( - Ytemp, Mathematics.svMult(0.5, K[0])))); - K[2] = Mathematics.svMult(h, DES.getValue(x + h / 2, Mathematics.vvAdd( - Ytemp, Mathematics.svMult(0.5, K[1])))); - K[3] = Mathematics.svMult(h, DES.getValue(x + h, Mathematics.vvAdd( - Ytemp, K[2]))); - - double[] change = Mathematics.svDiv(6, Mathematics.vvAdd(K[0], - Mathematics.vvAdd(Mathematics.svMult(2, K[1]), Mathematics - .vvAdd(Mathematics.svMult(2, K[2]), K[3])))); - for (int k = 0; k < change.length; k++) { - if (Double.isNaN(change[k])) { - unstableFlag = true; - change[k] = 0; - // return result; - } - } - return change; - } - - /** - * Linearized code for speed-up (no allocations). - * - * @param DES - * @param h - * @param x - * @param Ytemp - * @return - */ - public void rkTerm2(DESystem DES, double h, double x, double[] Ytemp, - double[] res) { - if (kVals == null) { // "static" vectors which are allocated only once - k0tmp = new double[DES.getDESystemDimension()]; - k1tmp = new double[DES.getDESystemDimension()]; - k2tmp = new double[DES.getDESystemDimension()]; - kVals = new double[4][DES.getDESystemDimension()]; - } - - // double[][] K = new double[4][]; - DES.getValue(x, Ytemp, kVals[0]); - Mathematics.svMult(h, kVals[0], kVals[0]); - - // K[0] = svMult(h, DES.getValue(x, Ytemp)); - - Mathematics.svMult(0.5, kVals[0], k0tmp); - Mathematics.vvAdd(Ytemp, k0tmp, k0tmp); - DES.getValue(x + h / 2, k0tmp, kVals[1]); - Mathematics.svMult(h, kVals[1], kVals[1]); - - // K[1] = svMult(h, DES.getValue(x + h / 2, vvAdd(Ytemp, svMult(0.5, - // K[0])))); - - Mathematics.svMult(0.5, kVals[1], k1tmp); - Mathematics.vvAdd(Ytemp, k1tmp, k1tmp); - DES.getValue(x + h / 2, k1tmp, kVals[2]); - Mathematics.svMult(h, kVals[2], kVals[2]); - - // K[2] = svMult(h, DES.getValue(x + h / 2, vvAdd(Ytemp, svMult(0.5, - // K[1])))); - - Mathematics.vvAdd(Ytemp, kVals[2], k2tmp); - DES.getValue(x + h, k2tmp, k1tmp); - Mathematics.svMult(h, k1tmp, kVals[3]); - - // K[3] = svMult(h, DES.getValue(x + h, vvAdd(Ytemp, K[2]))); - - Mathematics.svMult(2, kVals[2], k0tmp); - Mathematics.vvAdd(k0tmp, kVals[3], k0tmp); - - Mathematics.svMult(2, kVals[1], k1tmp); - Mathematics.vvAdd(k1tmp, k0tmp, k2tmp); - - Mathematics.vvAdd(kVals[0], k2tmp, k1tmp); - Mathematics.svDiv(6, k1tmp, res); - - // double[] change = svDiv(6, vvAdd(K[0], vvAdd(svMult(2, K[1]), - // vvAdd(svMult(2, K[2]), K[3])))); - // for (int i=0; i 0.00000001) System.out.println("!!! "); - // } - - // double[] change = svdiv(6, vvadd(kVals[0], vvadd(svmult(2, kVals[1]), - // vvadd(svmult(2, kVals[2]), kVals[3])))); - for (int k = 0; k < res.length; k++) { - if (Double.isNaN(res[k])) { - unstableFlag = true; - res[k] = 0; - // return result; - } - } - } - /** * @param stepSize */ @@ -254,9 +152,10 @@ public class RKSolver implements DESSolver, Serializable { return solveAtTimePoints(DES, initialValues, timePoints, true); } - /** - * - */ + /* + * (non-Javadoc) + * @see eva2.tools.math.des.DESSolver#solveAtTimePointsWithInitialConditions(eva2.tools.math.des.DESystem, double[][], double[]) + */ public double[][] solveAtTimePointsWithInitialConditions(DESystem DES, double[][] initConditions, double[] timePoints) { int order = DES.getDESystemDimension(); @@ -319,6 +218,109 @@ public class RKSolver implements DESSolver, Serializable { return solveByStepSize(DES, initialValues, timeBegin, timeEnd, true); } + /** + * @param DES + * @param h + * @param x + * @param Ytemp + * @return + */ + private double[] rkTerm(DESystem DES, double h, double x, double[] Ytemp) { + double[][] K = new double[4][]; + K[0] = Mathematics.svMult(h, DES.getValue(x, Ytemp)); + K[1] = Mathematics.svMult(h, DES.getValue(x + h / 2, Mathematics.vvAdd( + Ytemp, Mathematics.svMult(0.5, K[0])))); + K[2] = Mathematics.svMult(h, DES.getValue(x + h / 2, Mathematics.vvAdd( + Ytemp, Mathematics.svMult(0.5, K[1])))); + K[3] = Mathematics.svMult(h, DES.getValue(x + h, Mathematics.vvAdd( + Ytemp, K[2]))); + + double[] change = Mathematics.svDiv(6, Mathematics.vvAdd(K[0], + Mathematics.vvAdd(Mathematics.svMult(2, K[1]), Mathematics + .vvAdd(Mathematics.svMult(2, K[2]), K[3])))); + for (int k = 0; k < change.length; k++) { + if (Double.isNaN(change[k])) { + unstableFlag = true; + change[k] = 0; + // return result; + } + } + return change; + } + + /** + * Linearized code for speed-up (no allocations). + * + * @param DES + * @param h + * @param x + * @param Ytemp + * @return + */ + private void rkTerm2(DESystem DES, double h, double x, double[] Ytemp, + double[] res) { + if (kVals == null) { // "static" vectors which are allocated only once + k0tmp = new double[DES.getDESystemDimension()]; + k1tmp = new double[DES.getDESystemDimension()]; + k2tmp = new double[DES.getDESystemDimension()]; + kVals = new double[4][DES.getDESystemDimension()]; + } + + // double[][] K = new double[4][]; + DES.getValue(x, Ytemp, kVals[0]); + Mathematics.svMult(h, kVals[0], kVals[0]); + + // K[0] = svMult(h, DES.getValue(x, Ytemp)); + + Mathematics.svMult(0.5, kVals[0], k0tmp); + Mathematics.vvAdd(Ytemp, k0tmp, k0tmp); + DES.getValue(x + h / 2, k0tmp, kVals[1]); + Mathematics.svMult(h, kVals[1], kVals[1]); + + // K[1] = svMult(h, DES.getValue(x + h / 2, vvAdd(Ytemp, svMult(0.5, + // K[0])))); + + Mathematics.svMult(0.5, kVals[1], k1tmp); + Mathematics.vvAdd(Ytemp, k1tmp, k1tmp); + DES.getValue(x + h / 2, k1tmp, kVals[2]); + Mathematics.svMult(h, kVals[2], kVals[2]); + + // K[2] = svMult(h, DES.getValue(x + h / 2, vvAdd(Ytemp, svMult(0.5, + // K[1])))); + + Mathematics.vvAdd(Ytemp, kVals[2], k2tmp); + DES.getValue(x + h, k2tmp, k1tmp); + Mathematics.svMult(h, k1tmp, kVals[3]); + + // K[3] = svMult(h, DES.getValue(x + h, vvAdd(Ytemp, K[2]))); + + Mathematics.svMult(2, kVals[2], k0tmp); + Mathematics.vvAdd(k0tmp, kVals[3], k0tmp); + + Mathematics.svMult(2, kVals[1], k1tmp); + Mathematics.vvAdd(k1tmp, k0tmp, k2tmp); + + Mathematics.vvAdd(kVals[0], k2tmp, k1tmp); + Mathematics.svDiv(6, k1tmp, res); + + // double[] change = svDiv(6, vvAdd(K[0], vvAdd(svMult(2, K[1]), + // vvAdd(svMult(2, K[2]), K[3])))); + // for (int i=0; i 0.00000001) System.out.println("!!! "); + // } + + // double[] change = svdiv(6, vvadd(kVals[0], vvadd(svmult(2, kVals[1]), + // vvadd(svmult(2, kVals[2]), kVals[3])))); + for (int k = 0; k < res.length; k++) { + if (Double.isNaN(res[k])) { + unstableFlag = true; + res[k] = 0; + // return result; + } + } + } + /** * When set to TRUE, includeTimes will make the * solver to return a matrix with the first column containing the times. By