wrap lines *after* operator, to keep checkstyle happy
git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@613600 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
parent
ab1f56eacb
commit
540babda6c
|
@ -124,25 +124,25 @@ public class MessagesResources_fr
|
||||||
// org.apache.commons.math.ode.AdaptiveStepsizeIntegrator
|
// org.apache.commons.math.ode.AdaptiveStepsizeIntegrator
|
||||||
{ "minimal step size ({0}) reached, integration needs {1}",
|
{ "minimal step size ({0}) reached, integration needs {1}",
|
||||||
"pas minimal ({0}) atteint, l''int\u00e9gration n\u00e9cessite {1}" },
|
"pas minimal ({0}) atteint, l''int\u00e9gration n\u00e9cessite {1}" },
|
||||||
{ "dimensions mismatch: state vector has dimension {0},"
|
{ "dimensions mismatch: state vector has dimension {0}," +
|
||||||
+ " absolute tolerance vector has dimension {1}",
|
" absolute tolerance vector has dimension {1}",
|
||||||
"incompatibilit\u00e9 de dimensions entre le vecteur d''\u00e9tat ({0}),"
|
"incompatibilit\u00e9 de dimensions entre le vecteur d''\u00e9tat ({0})," +
|
||||||
+ " et le vecteur de tol\u00e9rance absolue ({1})" },
|
" et le vecteur de tol\u00e9rance absolue ({1})" },
|
||||||
{ "dimensions mismatch: state vector has dimension {0},"
|
{ "dimensions mismatch: state vector has dimension {0}," +
|
||||||
+ " relative tolerance vector has dimension {1}",
|
" relative tolerance vector has dimension {1}",
|
||||||
"incompatibilit\u00e9 de dimensions entre le vecteur d''\u00e9tat ({0}),"
|
"incompatibilit\u00e9 de dimensions entre le vecteur d''\u00e9tat ({0})," +
|
||||||
+ " et le vecteur de tol\u00e9rance relative ({1})" },
|
" et le vecteur de tol\u00e9rance relative ({1})" },
|
||||||
|
|
||||||
// org.apache.commons.math.ode.AdaptiveStepsizeIntegrator,
|
// org.apache.commons.math.ode.AdaptiveStepsizeIntegrator,
|
||||||
// org.apache.commons.math.ode.RungeKuttaIntegrator
|
// org.apache.commons.math.ode.RungeKuttaIntegrator
|
||||||
{ "dimensions mismatch: ODE problem has dimension {0},"
|
{ "dimensions mismatch: ODE problem has dimension {0}," +
|
||||||
+ " initial state vector has dimension {1}",
|
" initial state vector has dimension {1}",
|
||||||
"incompatibilit\u00e9 de dimensions entre le probl\u00e8me ODE ({0}),"
|
"incompatibilit\u00e9 de dimensions entre le probl\u00e8me ODE ({0})," +
|
||||||
+ " et le vecteur d''\u00e9tat initial ({1})" },
|
" et le vecteur d''\u00e9tat initial ({1})" },
|
||||||
{ "dimensions mismatch: ODE problem has dimension {0},"
|
{ "dimensions mismatch: ODE problem has dimension {0}," +
|
||||||
+ " final state vector has dimension {1}",
|
" final state vector has dimension {1}",
|
||||||
"incompatibilit\u00e9 de dimensions entre le probl\u00e8me ODE ({0}),"
|
"incompatibilit\u00e9 de dimensions entre le probl\u00e8me ODE ({0})," +
|
||||||
+ " et le vecteur d''\u00e9tat final ({1})" },
|
" et le vecteur d''\u00e9tat final ({1})" },
|
||||||
{ "too small integration interval: length = {0}",
|
{ "too small integration interval: length = {0}",
|
||||||
"intervalle d''int\u00e9gration trop petit : {0}" },
|
"intervalle d''int\u00e9gration trop petit : {0}" },
|
||||||
|
|
||||||
|
|
|
@ -65,9 +65,9 @@ public class BrentSolver extends UnivariateRealSolverImpl {
|
||||||
throws MaxIterationsExceededException, FunctionEvaluationException {
|
throws MaxIterationsExceededException, FunctionEvaluationException {
|
||||||
|
|
||||||
if (((initial - min) * (max -initial)) < 0) {
|
if (((initial - min) * (max -initial)) < 0) {
|
||||||
throw new IllegalArgumentException("Initial guess is not in search"
|
throw new IllegalArgumentException("Initial guess is not in search" +
|
||||||
+ " interval." + " Initial: " + initial
|
" interval." + " Initial: " + initial +
|
||||||
+ " Endpoints: [" + min + "," + max + "]");
|
" Endpoints: [" + min + "," + max + "]");
|
||||||
}
|
}
|
||||||
|
|
||||||
// return the initial guess if it is good enough
|
// return the initial guess if it is good enough
|
||||||
|
|
|
@ -146,10 +146,10 @@ public class PascalDistributionImpl extends AbstractIntegerDistribution
|
||||||
if (x < 0) {
|
if (x < 0) {
|
||||||
ret = 0.0;
|
ret = 0.0;
|
||||||
} else {
|
} else {
|
||||||
ret = MathUtils.binomialCoefficientDouble(x
|
ret = MathUtils.binomialCoefficientDouble(x +
|
||||||
+ getNumberOfSuccesses() - 1, getNumberOfSuccesses() - 1)
|
getNumberOfSuccesses() - 1, getNumberOfSuccesses() - 1) *
|
||||||
* Math.pow(getProbabilityOfSuccess(), getNumberOfSuccesses())
|
Math.pow(getProbabilityOfSuccess(), getNumberOfSuccesses()) *
|
||||||
* Math.pow(1.0 - getProbabilityOfSuccess(), x);
|
Math.pow(1.0 - getProbabilityOfSuccess(), x);
|
||||||
}
|
}
|
||||||
return ret;
|
return ret;
|
||||||
}
|
}
|
||||||
|
|
|
@ -164,9 +164,9 @@ public class GaussNewtonEstimator extends AbstractEstimator implements Serializa
|
||||||
previous = cost;
|
previous = cost;
|
||||||
updateResidualsAndCost();
|
updateResidualsAndCost();
|
||||||
|
|
||||||
} while ((getCostEvaluations() < 2)
|
} while ((getCostEvaluations() < 2) ||
|
||||||
|| (Math.abs(previous - cost) > (cost * steadyStateThreshold)
|
(Math.abs(previous - cost) > (cost * steadyStateThreshold) &&
|
||||||
&& (Math.abs(cost) > convergence)));
|
(Math.abs(cost) > convergence)));
|
||||||
|
|
||||||
}
|
}
|
||||||
|
|
||||||
|
|
|
@ -257,8 +257,7 @@ public class LevenbergMarquardtEstimator extends AbstractEstimator implements Se
|
||||||
xNorm = Math.sqrt(xNorm);
|
xNorm = Math.sqrt(xNorm);
|
||||||
|
|
||||||
// initialize the step bound delta
|
// initialize the step bound delta
|
||||||
delta = (xNorm == 0)
|
delta = (xNorm == 0) ? initialStepBoundFactor : (initialStepBoundFactor * xNorm);
|
||||||
? initialStepBoundFactor : (initialStepBoundFactor * xNorm);
|
|
||||||
|
|
||||||
}
|
}
|
||||||
|
|
||||||
|
@ -388,30 +387,28 @@ public class LevenbergMarquardtEstimator extends AbstractEstimator implements Se
|
||||||
}
|
}
|
||||||
|
|
||||||
// tests for convergence.
|
// tests for convergence.
|
||||||
if (((Math.abs(actRed) <= costRelativeTolerance)
|
if (((Math.abs(actRed) <= costRelativeTolerance) &&
|
||||||
&& (preRed <= costRelativeTolerance)
|
(preRed <= costRelativeTolerance) &&
|
||||||
&& (ratio <= 2.0))
|
(ratio <= 2.0)) ||
|
||||||
|| (delta <= parRelativeTolerance * xNorm)) {
|
(delta <= parRelativeTolerance * xNorm)) {
|
||||||
return;
|
return;
|
||||||
}
|
}
|
||||||
|
|
||||||
// tests for termination and stringent tolerances
|
// tests for termination and stringent tolerances
|
||||||
// (2.2204e-16 is the machine epsilon for IEEE754)
|
// (2.2204e-16 is the machine epsilon for IEEE754)
|
||||||
if ((Math.abs(actRed) <= 2.2204e-16)
|
if ((Math.abs(actRed) <= 2.2204e-16) && (preRed <= 2.2204e-16) && (ratio <= 2.0)) {
|
||||||
&& (preRed <= 2.2204e-16)
|
throw new EstimationException("cost relative tolerance is too small ({0})," +
|
||||||
&& (ratio <= 2.0)) {
|
" no further reduction in the" +
|
||||||
throw new EstimationException("cost relative tolerance is too small ({0}),"
|
" sum of squares is possible",
|
||||||
+ " no further reduction in the"
|
|
||||||
+ " sum of squares is possible",
|
|
||||||
new Object[] { new Double(costRelativeTolerance) });
|
new Object[] { new Double(costRelativeTolerance) });
|
||||||
} else if (delta <= 2.2204e-16 * xNorm) {
|
} else if (delta <= 2.2204e-16 * xNorm) {
|
||||||
throw new EstimationException("parameters relative tolerance is too small"
|
throw new EstimationException("parameters relative tolerance is too small" +
|
||||||
+ " ({0}), no further improvement in"
|
" ({0}), no further improvement in" +
|
||||||
+ " the approximate solution is possible",
|
" the approximate solution is possible",
|
||||||
new Object[] { new Double(parRelativeTolerance) });
|
new Object[] { new Double(parRelativeTolerance) });
|
||||||
} else if (maxCosine <= 2.2204e-16) {
|
} else if (maxCosine <= 2.2204e-16) {
|
||||||
throw new EstimationException("orthogonality tolerance is too small ({0}),"
|
throw new EstimationException("orthogonality tolerance is too small ({0})," +
|
||||||
+ " solution is orthogonal to the jacobian",
|
" solution is orthogonal to the jacobian",
|
||||||
new Object[] { new Double(orthoTolerance) });
|
new Object[] { new Double(orthoTolerance) });
|
||||||
}
|
}
|
||||||
|
|
||||||
|
@ -553,8 +550,8 @@ public class LevenbergMarquardtEstimator extends AbstractEstimator implements Se
|
||||||
|
|
||||||
// if the function is small enough, accept the current value
|
// if the function is small enough, accept the current value
|
||||||
// of lmPar, also test for the exceptional cases where parl is zero
|
// of lmPar, also test for the exceptional cases where parl is zero
|
||||||
if ((Math.abs(fp) <= 0.1 * delta)
|
if ((Math.abs(fp) <= 0.1 * delta) ||
|
||||||
|| ((parl == 0) && (fp <= previousFP) && (previousFP < 0))) {
|
((parl == 0) && (fp <= previousFP) && (previousFP < 0))) {
|
||||||
return;
|
return;
|
||||||
}
|
}
|
||||||
|
|
||||||
|
|
|
@ -192,10 +192,10 @@ public class Rotation implements Serializable {
|
||||||
throws NotARotationMatrixException {
|
throws NotARotationMatrixException {
|
||||||
|
|
||||||
// dimension check
|
// dimension check
|
||||||
if ((m.length != 3) || (m[0].length != 3)
|
if ((m.length != 3) || (m[0].length != 3) ||
|
||||||
|| (m[1].length != 3) || (m[2].length != 3)) {
|
(m[1].length != 3) || (m[2].length != 3)) {
|
||||||
throw new NotARotationMatrixException("a {0}x{1} matrix"
|
throw new NotARotationMatrixException("a {0}x{1} matrix" +
|
||||||
+ " cannot be a rotation matrix",
|
" cannot be a rotation matrix",
|
||||||
new String[] {
|
new String[] {
|
||||||
Integer.toString(m.length),
|
Integer.toString(m.length),
|
||||||
Integer.toString(m[0].length)
|
Integer.toString(m[0].length)
|
||||||
|
@ -206,12 +206,12 @@ public class Rotation implements Serializable {
|
||||||
double[][] ort = orthogonalizeMatrix(m, threshold);
|
double[][] ort = orthogonalizeMatrix(m, threshold);
|
||||||
|
|
||||||
// check the sign of the determinant
|
// check the sign of the determinant
|
||||||
double det = ort[0][0] * (ort[1][1] * ort[2][2] - ort[2][1] * ort[1][2])
|
double det = ort[0][0] * (ort[1][1] * ort[2][2] - ort[2][1] * ort[1][2]) -
|
||||||
- ort[1][0] * (ort[0][1] * ort[2][2] - ort[2][1] * ort[0][2])
|
ort[1][0] * (ort[0][1] * ort[2][2] - ort[2][1] * ort[0][2]) +
|
||||||
+ ort[2][0] * (ort[0][1] * ort[1][2] - ort[1][1] * ort[0][2]);
|
ort[2][0] * (ort[0][1] * ort[1][2] - ort[1][1] * ort[0][2]);
|
||||||
if (det < 0.0) {
|
if (det < 0.0) {
|
||||||
throw new NotARotationMatrixException("the closest orthogonal matrix"
|
throw new NotARotationMatrixException("the closest orthogonal matrix" +
|
||||||
+ " has a negative determinant {0}",
|
" has a negative determinant {0}",
|
||||||
new String[] {
|
new String[] {
|
||||||
Double.toString(det)
|
Double.toString(det)
|
||||||
});
|
});
|
||||||
|
@ -337,9 +337,9 @@ public class Rotation implements Serializable {
|
||||||
Vector3D k = new Vector3D(dy1 * dz2 - dz1 * dy2,
|
Vector3D k = new Vector3D(dy1 * dz2 - dz1 * dy2,
|
||||||
dz1 * dx2 - dx1 * dz2,
|
dz1 * dx2 - dx1 * dz2,
|
||||||
dx1 * dy2 - dy1 * dx2);
|
dx1 * dy2 - dy1 * dx2);
|
||||||
double c = k.getX() * (u1y * u2z - u1z * u2y)
|
double c = k.getX() * (u1y * u2z - u1z * u2y) +
|
||||||
+ k.getY() * (u1z * u2x - u1x * u2z)
|
k.getY() * (u1z * u2x - u1x * u2z) +
|
||||||
+ k.getZ() * (u1x * u2y - u1y * u2x);
|
k.getZ() * (u1x * u2y - u1y * u2x);
|
||||||
|
|
||||||
if (c == 0) {
|
if (c == 0) {
|
||||||
// the (q1, q2, q3) vector is in the (u1, u2) plane
|
// the (q1, q2, q3) vector is in the (u1, u2) plane
|
||||||
|
@ -359,9 +359,9 @@ public class Rotation implements Serializable {
|
||||||
k = new Vector3D(dy1 * dz3 - dz1 * dy3,
|
k = new Vector3D(dy1 * dz3 - dz1 * dy3,
|
||||||
dz1 * dx3 - dx1 * dz3,
|
dz1 * dx3 - dx1 * dz3,
|
||||||
dx1 * dy3 - dy1 * dx3);
|
dx1 * dy3 - dy1 * dx3);
|
||||||
c = k.getX() * (u1y * u3z - u1z * u3y)
|
c = k.getX() * (u1y * u3z - u1z * u3y) +
|
||||||
+ k.getY() * (u1z * u3x - u1x * u3z)
|
k.getY() * (u1z * u3x - u1x * u3z) +
|
||||||
+ k.getZ() * (u1x * u3y - u1y * u3x);
|
k.getZ() * (u1x * u3y - u1y * u3x);
|
||||||
|
|
||||||
if (c == 0) {
|
if (c == 0) {
|
||||||
// the (q1, q2, q3) vector is aligned with u1:
|
// the (q1, q2, q3) vector is aligned with u1:
|
||||||
|
@ -369,9 +369,9 @@ public class Rotation implements Serializable {
|
||||||
k = new Vector3D(dy2 * dz3 - dz2 * dy3,
|
k = new Vector3D(dy2 * dz3 - dz2 * dy3,
|
||||||
dz2 * dx3 - dx2 * dz3,
|
dz2 * dx3 - dx2 * dz3,
|
||||||
dx2 * dy3 - dy2 * dx3);
|
dx2 * dy3 - dy2 * dx3);
|
||||||
c = k.getX() * (u2y * u3z - u2z * u3y)
|
c = k.getX() * (u2y * u3z - u2z * u3y) +
|
||||||
+ k.getY() * (u2z * u3x - u2x * u3z)
|
k.getY() * (u2z * u3x - u2x * u3z) +
|
||||||
+ k.getZ() * (u2x * u3y - u2y * u3x);
|
k.getZ() * (u2x * u3y - u2y * u3x);
|
||||||
|
|
||||||
if (c == 0) {
|
if (c == 0) {
|
||||||
// the (q1, q2, q3) vector is aligned with everything
|
// the (q1, q2, q3) vector is aligned with everything
|
||||||
|
@ -986,9 +986,9 @@ public class Rotation implements Serializable {
|
||||||
double corr22 = o2[2] - m2[2];
|
double corr22 = o2[2] - m2[2];
|
||||||
|
|
||||||
// Frobenius norm of the correction
|
// Frobenius norm of the correction
|
||||||
fn1 = corr00 * corr00 + corr01 * corr01 + corr02 * corr02
|
fn1 = corr00 * corr00 + corr01 * corr01 + corr02 * corr02 +
|
||||||
+ corr10 * corr10 + corr11 * corr11 + corr12 * corr12
|
corr10 * corr10 + corr11 * corr11 + corr12 * corr12 +
|
||||||
+ corr20 * corr20 + corr21 * corr21 + corr22 * corr22;
|
corr20 * corr20 + corr21 * corr21 + corr22 * corr22;
|
||||||
|
|
||||||
// convergence test
|
// convergence test
|
||||||
if (Math.abs(fn1 - fn) <= threshold)
|
if (Math.abs(fn1 - fn) <= threshold)
|
||||||
|
@ -1009,8 +1009,8 @@ public class Rotation implements Serializable {
|
||||||
}
|
}
|
||||||
|
|
||||||
// the algorithm did not converge after 10 iterations
|
// the algorithm did not converge after 10 iterations
|
||||||
throw new NotARotationMatrixException("unable to orthogonalize matrix"
|
throw new NotARotationMatrixException("unable to orthogonalize matrix" +
|
||||||
+ " in {0} iterations",
|
" in {0} iterations",
|
||||||
new String[] {
|
new String[] {
|
||||||
Integer.toString(i - 1)
|
Integer.toString(i - 1)
|
||||||
});
|
});
|
||||||
|
|
|
@ -95,10 +95,10 @@ class ClassicalRungeKuttaStepInterpolator
|
||||||
double coeff4 = s * ((-fourTheta - 1) * theta - 1);
|
double coeff4 = s * ((-fourTheta - 1) * theta - 1);
|
||||||
|
|
||||||
for (int i = 0; i < interpolatedState.length; ++i) {
|
for (int i = 0; i < interpolatedState.length; ++i) {
|
||||||
interpolatedState[i] = currentState[i]
|
interpolatedState[i] = currentState[i] +
|
||||||
+ coeff1 * yDotK[0][i]
|
coeff1 * yDotK[0][i] +
|
||||||
+ coeff23 * (yDotK[1][i] + yDotK[2][i])
|
coeff23 * (yDotK[1][i] + yDotK[2][i]) +
|
||||||
+ coeff4 * yDotK[3][i];
|
coeff4 * yDotK[3][i];
|
||||||
}
|
}
|
||||||
|
|
||||||
}
|
}
|
||||||
|
|
|
@ -286,10 +286,10 @@ public class ContinuousOutputModel
|
||||||
double dt1 = time - tMax;
|
double dt1 = time - tMax;
|
||||||
double dt2 = time - tMed;
|
double dt2 = time - tMed;
|
||||||
double dt3 = time - tMin;
|
double dt3 = time - tMin;
|
||||||
double iLagrange = ( (dt2 * dt3 * d23) * iMax
|
double iLagrange = ((dt2 * dt3 * d23) * iMax -
|
||||||
- (dt1 * dt3 * d13) * iMed
|
(dt1 * dt3 * d13) * iMed +
|
||||||
+ (dt1 * dt2 * d12) * iMin)
|
(dt1 * dt2 * d12) * iMin) /
|
||||||
/ (d12 * d23 * d13);
|
(d12 * d23 * d13);
|
||||||
index = (int) Math.rint(iLagrange);
|
index = (int) Math.rint(iLagrange);
|
||||||
}
|
}
|
||||||
|
|
||||||
|
@ -306,8 +306,8 @@ public class ContinuousOutputModel
|
||||||
|
|
||||||
// now the table slice is very small, we perform an iterative search
|
// now the table slice is very small, we perform an iterative search
|
||||||
index = iMin;
|
index = iMin;
|
||||||
while ((index <= iMax)
|
while ((index <= iMax) &&
|
||||||
&& (locatePoint(time, (StepInterpolator) steps.get(index)) > 0)) {
|
(locatePoint(time, (StepInterpolator) steps.get(index)) > 0)) {
|
||||||
++index;
|
++index;
|
||||||
}
|
}
|
||||||
|
|
||||||
|
@ -316,8 +316,8 @@ public class ContinuousOutputModel
|
||||||
si.setInterpolatedTime(time);
|
si.setInterpolatedTime(time);
|
||||||
|
|
||||||
} catch (DerivativeException de) {
|
} catch (DerivativeException de) {
|
||||||
throw new RuntimeException("unexpected DerivativeException caught: "
|
throw new RuntimeException("unexpected DerivativeException caught: " +
|
||||||
+ de.getMessage());
|
de.getMessage());
|
||||||
}
|
}
|
||||||
|
|
||||||
}
|
}
|
||||||
|
|
|
@ -132,14 +132,14 @@ public class DormandPrince54Integrator
|
||||||
double error = 0;
|
double error = 0;
|
||||||
|
|
||||||
for (int j = 0; j < y0.length; ++j) {
|
for (int j = 0; j < y0.length; ++j) {
|
||||||
double errSum = e1 * yDotK[0][j] + e3 * yDotK[2][j]
|
double errSum = e1 * yDotK[0][j] + e3 * yDotK[2][j] +
|
||||||
+ e4 * yDotK[3][j] + e5 * yDotK[4][j]
|
e4 * yDotK[3][j] + e5 * yDotK[4][j] +
|
||||||
+ e6 * yDotK[5][j] + e7 * yDotK[6][j];
|
e6 * yDotK[5][j] + e7 * yDotK[6][j];
|
||||||
|
|
||||||
double yScale = Math.max(Math.abs(y0[j]), Math.abs(y1[j]));
|
double yScale = Math.max(Math.abs(y0[j]), Math.abs(y1[j]));
|
||||||
double tol = (vecAbsoluteTolerance == null)
|
double tol = (vecAbsoluteTolerance == null) ?
|
||||||
? (scalAbsoluteTolerance + scalRelativeTolerance * yScale)
|
(scalAbsoluteTolerance + scalRelativeTolerance * yScale) :
|
||||||
: (vecAbsoluteTolerance[j] + vecRelativeTolerance[j] * yScale);
|
(vecAbsoluteTolerance[j] + vecRelativeTolerance[j] * yScale);
|
||||||
double ratio = h * errSum / tol;
|
double ratio = h * errSum / tol;
|
||||||
error += ratio * ratio;
|
error += ratio * ratio;
|
||||||
|
|
||||||
|
|
|
@ -135,12 +135,12 @@ class DormandPrince54StepInterpolator
|
||||||
|
|
||||||
// we need to compute the interpolation vectors for this time step
|
// we need to compute the interpolation vectors for this time step
|
||||||
for (int i = 0; i < interpolatedState.length; ++i) {
|
for (int i = 0; i < interpolatedState.length; ++i) {
|
||||||
v1[i] = h * (a70 * yDotK[0][i] + a72 * yDotK[2][i] + a73 * yDotK[3][i]
|
v1[i] = h * (a70 * yDotK[0][i] + a72 * yDotK[2][i] + a73 * yDotK[3][i] +
|
||||||
+ a74 * yDotK[4][i] + a75 * yDotK[5][i]);
|
a74 * yDotK[4][i] + a75 * yDotK[5][i]);
|
||||||
v2[i] = h * yDotK[0][i] - v1[i];
|
v2[i] = h * yDotK[0][i] - v1[i];
|
||||||
v3[i] = v1[i] - v2[i] - h * yDotK[6][i];
|
v3[i] = v1[i] - v2[i] - h * yDotK[6][i];
|
||||||
v4[i] = h * (d0 * yDotK[0][i] + d2 * yDotK[2][i] + d3 * yDotK[3][i]
|
v4[i] = h * (d0 * yDotK[0][i] + d2 * yDotK[2][i] + d3 * yDotK[3][i] +
|
||||||
+ d4 * yDotK[4][i] + d5 * yDotK[5][i] + d6 * yDotK[6][i]);
|
d4 * yDotK[4][i] + d5 * yDotK[5][i] + d6 * yDotK[6][i]);
|
||||||
}
|
}
|
||||||
|
|
||||||
vectorsInitialized = true;
|
vectorsInitialized = true;
|
||||||
|
@ -150,11 +150,8 @@ class DormandPrince54StepInterpolator
|
||||||
// interpolate
|
// interpolate
|
||||||
double eta = oneMinusThetaH / h;
|
double eta = oneMinusThetaH / h;
|
||||||
for (int i = 0; i < interpolatedState.length; ++i) {
|
for (int i = 0; i < interpolatedState.length; ++i) {
|
||||||
interpolatedState[i] = currentState[i]
|
interpolatedState[i] = currentState[i] -
|
||||||
- eta * (v1[i]
|
eta * (v1[i] - theta * (v2[i] + theta * (v3[i] + eta * v4[i])));
|
||||||
- theta * (v2[i]
|
|
||||||
+ theta * (v3[i]
|
|
||||||
+ eta * v4[i])));
|
|
||||||
}
|
}
|
||||||
|
|
||||||
}
|
}
|
||||||
|
|
|
@ -226,19 +226,19 @@ public class DormandPrince853Integrator
|
||||||
double error2 = 0;
|
double error2 = 0;
|
||||||
|
|
||||||
for (int j = 0; j < y0.length; ++j) {
|
for (int j = 0; j < y0.length; ++j) {
|
||||||
double errSum1 = e1_01 * yDotK[0][j] + e1_06 * yDotK[5][j]
|
double errSum1 = e1_01 * yDotK[0][j] + e1_06 * yDotK[5][j] +
|
||||||
+ e1_07 * yDotK[6][j] + e1_08 * yDotK[7][j]
|
e1_07 * yDotK[6][j] + e1_08 * yDotK[7][j] +
|
||||||
+ e1_09 * yDotK[8][j] + e1_10 * yDotK[9][j]
|
e1_09 * yDotK[8][j] + e1_10 * yDotK[9][j] +
|
||||||
+ e1_11 * yDotK[10][j] + e1_12 * yDotK[11][j];
|
e1_11 * yDotK[10][j] + e1_12 * yDotK[11][j];
|
||||||
double errSum2 = e2_01 * yDotK[0][j] + e2_06 * yDotK[5][j]
|
double errSum2 = e2_01 * yDotK[0][j] + e2_06 * yDotK[5][j] +
|
||||||
+ e2_07 * yDotK[6][j] + e2_08 * yDotK[7][j]
|
e2_07 * yDotK[6][j] + e2_08 * yDotK[7][j] +
|
||||||
+ e2_09 * yDotK[8][j] + e2_10 * yDotK[9][j]
|
e2_09 * yDotK[8][j] + e2_10 * yDotK[9][j] +
|
||||||
+ e2_11 * yDotK[10][j] + e2_12 * yDotK[11][j];
|
e2_11 * yDotK[10][j] + e2_12 * yDotK[11][j];
|
||||||
|
|
||||||
double yScale = Math.max(Math.abs(y0[j]), Math.abs(y1[j]));
|
double yScale = Math.max(Math.abs(y0[j]), Math.abs(y1[j]));
|
||||||
double tol = (vecAbsoluteTolerance == null)
|
double tol = (vecAbsoluteTolerance == null) ?
|
||||||
? (scalAbsoluteTolerance + scalRelativeTolerance * yScale)
|
(scalAbsoluteTolerance + scalRelativeTolerance * yScale) :
|
||||||
: (vecAbsoluteTolerance[j] + vecRelativeTolerance[j] * yScale);
|
(vecAbsoluteTolerance[j] + vecRelativeTolerance[j] * yScale);
|
||||||
double ratio1 = errSum1 / tol;
|
double ratio1 = errSum1 / tol;
|
||||||
error1 += ratio1 * ratio1;
|
error1 += ratio1 * ratio1;
|
||||||
double ratio2 = errSum2 / tol;
|
double ratio2 = errSum2 / tol;
|
||||||
|
|
|
@ -172,18 +172,18 @@ class DormandPrince853StepInterpolator
|
||||||
|
|
||||||
// compute the interpolation vectors for this time step
|
// compute the interpolation vectors for this time step
|
||||||
for (int i = 0; i < interpolatedState.length; ++i) {
|
for (int i = 0; i < interpolatedState.length; ++i) {
|
||||||
v[0][i] = h * (b_01 * yDotK[0][i] + b_06 * yDotK[5][i] + b_07 * yDotK[6][i]
|
v[0][i] = h * (b_01 * yDotK[0][i] + b_06 * yDotK[5][i] + b_07 * yDotK[6][i] +
|
||||||
+ b_08 * yDotK[7][i] + b_09 * yDotK[8][i] + b_10 * yDotK[9][i]
|
b_08 * yDotK[7][i] + b_09 * yDotK[8][i] + b_10 * yDotK[9][i] +
|
||||||
+ b_11 * yDotK[10][i] + b_12 * yDotK[11][i]);
|
b_11 * yDotK[10][i] + b_12 * yDotK[11][i]);
|
||||||
v[1][i] = h * yDotK[0][i] - v[0][i];
|
v[1][i] = h * yDotK[0][i] - v[0][i];
|
||||||
v[2][i] = v[0][i] - v[1][i] - h * yDotK[12][i];
|
v[2][i] = v[0][i] - v[1][i] - h * yDotK[12][i];
|
||||||
for (int k = 0; k < d.length; ++k) {
|
for (int k = 0; k < d.length; ++k) {
|
||||||
v[k+3][i] = h * (d[k][0] * yDotK[0][i] + d[k][1] * yDotK[5][i] + d[k][2] * yDotK[6][i]
|
v[k+3][i] = h * (d[k][0] * yDotK[0][i] + d[k][1] * yDotK[5][i] + d[k][2] * yDotK[6][i] +
|
||||||
+ d[k][3] * yDotK[7][i] + d[k][4] * yDotK[8][i] + d[k][5] * yDotK[9][i]
|
d[k][3] * yDotK[7][i] + d[k][4] * yDotK[8][i] + d[k][5] * yDotK[9][i] +
|
||||||
+ d[k][6] * yDotK[10][i] + d[k][7] * yDotK[11][i] + d[k][8] * yDotK[12][i]
|
d[k][6] * yDotK[10][i] + d[k][7] * yDotK[11][i] + d[k][8] * yDotK[12][i] +
|
||||||
+ d[k][9] * yDotKLast[0][i]
|
d[k][9] * yDotKLast[0][i] +
|
||||||
+ d[k][10] * yDotKLast[1][i]
|
d[k][10] * yDotKLast[1][i] +
|
||||||
+ d[k][11] * yDotKLast[2][i]);
|
d[k][11] * yDotKLast[2][i]);
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
|
@ -194,14 +194,10 @@ class DormandPrince853StepInterpolator
|
||||||
double eta = oneMinusThetaH / h;
|
double eta = oneMinusThetaH / h;
|
||||||
|
|
||||||
for (int i = 0; i < interpolatedState.length; ++i) {
|
for (int i = 0; i < interpolatedState.length; ++i) {
|
||||||
interpolatedState[i] = currentState[i]
|
interpolatedState[i] =
|
||||||
- eta * (v[0][i]
|
currentState[i] - eta * (v[0][i] - theta * (v[1][i] +
|
||||||
- theta * (v[1][i]
|
theta * (v[2][i] + eta * (v[3][i] + theta * (v[4][i] +
|
||||||
+ theta * (v[2][i]
|
eta * (v[5][i] + theta * (v[6][i])))))));
|
||||||
+ eta * (v[3][i]
|
|
||||||
+ theta * (v[4][i]
|
|
||||||
+ eta * (v[5][i]
|
|
||||||
+ theta * (v[6][i])))))));
|
|
||||||
}
|
}
|
||||||
|
|
||||||
}
|
}
|
||||||
|
@ -225,29 +221,29 @@ class DormandPrince853StepInterpolator
|
||||||
|
|
||||||
// k14
|
// k14
|
||||||
for (int j = 0; j < currentState.length; ++j) {
|
for (int j = 0; j < currentState.length; ++j) {
|
||||||
s = k14_01 * yDotK[0][j] + k14_06 * yDotK[5][j] + k14_07 * yDotK[6][j]
|
s = k14_01 * yDotK[0][j] + k14_06 * yDotK[5][j] + k14_07 * yDotK[6][j] +
|
||||||
+ k14_08 * yDotK[7][j] + k14_09 * yDotK[8][j] + k14_10 * yDotK[9][j]
|
k14_08 * yDotK[7][j] + k14_09 * yDotK[8][j] + k14_10 * yDotK[9][j] +
|
||||||
+ k14_11 * yDotK[10][j] + k14_12 * yDotK[11][j] + k14_13 * yDotK[12][j];
|
k14_11 * yDotK[10][j] + k14_12 * yDotK[11][j] + k14_13 * yDotK[12][j];
|
||||||
yTmp[j] = currentState[j] + h * s;
|
yTmp[j] = currentState[j] + h * s;
|
||||||
}
|
}
|
||||||
equations.computeDerivatives(previousTime + c14 * h, yTmp, yDotKLast[0]);
|
equations.computeDerivatives(previousTime + c14 * h, yTmp, yDotKLast[0]);
|
||||||
|
|
||||||
// k15
|
// k15
|
||||||
for (int j = 0; j < currentState.length; ++j) {
|
for (int j = 0; j < currentState.length; ++j) {
|
||||||
s = k15_01 * yDotK[0][j] + k15_06 * yDotK[5][j] + k15_07 * yDotK[6][j]
|
s = k15_01 * yDotK[0][j] + k15_06 * yDotK[5][j] + k15_07 * yDotK[6][j] +
|
||||||
+ k15_08 * yDotK[7][j] + k15_09 * yDotK[8][j] + k15_10 * yDotK[9][j]
|
k15_08 * yDotK[7][j] + k15_09 * yDotK[8][j] + k15_10 * yDotK[9][j] +
|
||||||
+ k15_11 * yDotK[10][j] + k15_12 * yDotK[11][j] + k15_13 * yDotK[12][j]
|
k15_11 * yDotK[10][j] + k15_12 * yDotK[11][j] + k15_13 * yDotK[12][j] +
|
||||||
+ k15_14 * yDotKLast[0][j];
|
k15_14 * yDotKLast[0][j];
|
||||||
yTmp[j] = currentState[j] + h * s;
|
yTmp[j] = currentState[j] + h * s;
|
||||||
}
|
}
|
||||||
equations.computeDerivatives(previousTime + c15 * h, yTmp, yDotKLast[1]);
|
equations.computeDerivatives(previousTime + c15 * h, yTmp, yDotKLast[1]);
|
||||||
|
|
||||||
// k16
|
// k16
|
||||||
for (int j = 0; j < currentState.length; ++j) {
|
for (int j = 0; j < currentState.length; ++j) {
|
||||||
s = k16_01 * yDotK[0][j] + k16_06 * yDotK[5][j] + k16_07 * yDotK[6][j]
|
s = k16_01 * yDotK[0][j] + k16_06 * yDotK[5][j] + k16_07 * yDotK[6][j] +
|
||||||
+ k16_08 * yDotK[7][j] + k16_09 * yDotK[8][j] + k16_10 * yDotK[9][j]
|
k16_08 * yDotK[7][j] + k16_09 * yDotK[8][j] + k16_10 * yDotK[9][j] +
|
||||||
+ k16_11 * yDotK[10][j] + k16_12 * yDotK[11][j] + k16_13 * yDotK[12][j]
|
k16_11 * yDotK[10][j] + k16_12 * yDotK[11][j] + k16_13 * yDotK[12][j] +
|
||||||
+ k16_14 * yDotKLast[0][j] + k16_15 * yDotKLast[1][j];
|
k16_14 * yDotKLast[0][j] + k16_15 * yDotKLast[1][j];
|
||||||
yTmp[j] = currentState[j] + h * s;
|
yTmp[j] = currentState[j] + h * s;
|
||||||
}
|
}
|
||||||
equations.computeDerivatives(previousTime + c16 * h, yTmp, yDotKLast[2]);
|
equations.computeDerivatives(previousTime + c16 * h, yTmp, yDotKLast[2]);
|
||||||
|
|
|
@ -220,8 +220,8 @@ public abstract class EmbeddedRungeKuttaIntegrator
|
||||||
stepSize = hNew;
|
stepSize = hNew;
|
||||||
|
|
||||||
// step adjustment near bounds
|
// step adjustment near bounds
|
||||||
if ((forward && (stepStart + stepSize > t))
|
if ((forward && (stepStart + stepSize > t)) ||
|
||||||
|| ((! forward) && (stepStart + stepSize < t))) {
|
((! forward) && (stepStart + stepSize < t))) {
|
||||||
stepSize = t - stepStart;
|
stepSize = t - stepStart;
|
||||||
}
|
}
|
||||||
|
|
||||||
|
|
|
@ -97,9 +97,9 @@ class GillStepInterpolator
|
||||||
double coeff4 = s * (1 + theta * (1 + fourTheta));
|
double coeff4 = s * (1 + theta * (1 + fourTheta));
|
||||||
|
|
||||||
for (int i = 0; i < interpolatedState.length; ++i) {
|
for (int i = 0; i < interpolatedState.length; ++i) {
|
||||||
interpolatedState[i] = currentState[i]
|
interpolatedState[i] = currentState[i] -
|
||||||
- coeff1 * yDotK[0][i] - coeff2 * yDotK[1][i]
|
coeff1 * yDotK[0][i] - coeff2 * yDotK[1][i] -
|
||||||
- coeff3 * yDotK[2][i] - coeff4 * yDotK[3][i];
|
coeff3 * yDotK[2][i] - coeff4 * yDotK[3][i];
|
||||||
}
|
}
|
||||||
|
|
||||||
}
|
}
|
||||||
|
|
|
@ -104,8 +104,7 @@ public class GraggBulirschStoerIntegrator
|
||||||
double scalAbsoluteTolerance,
|
double scalAbsoluteTolerance,
|
||||||
double scalRelativeTolerance) {
|
double scalRelativeTolerance) {
|
||||||
super(minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance);
|
super(minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance);
|
||||||
denseOutput = (handler.requiresDenseOutput()
|
denseOutput = (handler.requiresDenseOutput() || (! switchesHandler.isEmpty()));
|
||||||
|| (! switchesHandler.isEmpty()));
|
|
||||||
setStabilityCheck(true, -1, -1, -1);
|
setStabilityCheck(true, -1, -1, -1);
|
||||||
setStepsizeControl(-1, -1, -1, -1);
|
setStepsizeControl(-1, -1, -1, -1);
|
||||||
setOrderControl(-1, -1, -1);
|
setOrderControl(-1, -1, -1);
|
||||||
|
@ -127,8 +126,7 @@ public class GraggBulirschStoerIntegrator
|
||||||
double[] vecAbsoluteTolerance,
|
double[] vecAbsoluteTolerance,
|
||||||
double[] vecRelativeTolerance) {
|
double[] vecRelativeTolerance) {
|
||||||
super(minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance);
|
super(minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance);
|
||||||
denseOutput = (handler.requiresDenseOutput()
|
denseOutput = (handler.requiresDenseOutput() || (! switchesHandler.isEmpty()));
|
||||||
|| (! switchesHandler.isEmpty()));
|
|
||||||
setStabilityCheck(true, -1, -1, -1);
|
setStabilityCheck(true, -1, -1, -1);
|
||||||
setStepsizeControl(-1, -1, -1, -1);
|
setStepsizeControl(-1, -1, -1, -1);
|
||||||
setOrderControl(-1, -1, -1);
|
setOrderControl(-1, -1, -1);
|
||||||
|
@ -277,8 +275,7 @@ public class GraggBulirschStoerIntegrator
|
||||||
public void setStepHandler (StepHandler handler) {
|
public void setStepHandler (StepHandler handler) {
|
||||||
|
|
||||||
super.setStepHandler(handler);
|
super.setStepHandler(handler);
|
||||||
denseOutput = (handler.requiresDenseOutput()
|
denseOutput = (handler.requiresDenseOutput() || (! switchesHandler.isEmpty()));
|
||||||
|| (! switchesHandler.isEmpty()));
|
|
||||||
|
|
||||||
// reinitialize the arrays
|
// reinitialize the arrays
|
||||||
initializeArrays();
|
initializeArrays();
|
||||||
|
@ -299,8 +296,7 @@ public class GraggBulirschStoerIntegrator
|
||||||
double convergence,
|
double convergence,
|
||||||
int maxIterationCount) {
|
int maxIterationCount) {
|
||||||
super.addSwitchingFunction(function, maxCheckInterval, convergence, maxIterationCount);
|
super.addSwitchingFunction(function, maxCheckInterval, convergence, maxIterationCount);
|
||||||
denseOutput = (handler.requiresDenseOutput()
|
denseOutput = (handler.requiresDenseOutput() || (! switchesHandler.isEmpty()));
|
||||||
|| (! switchesHandler.isEmpty()));
|
|
||||||
|
|
||||||
// reinitialize the arrays
|
// reinitialize the arrays
|
||||||
initializeArrays();
|
initializeArrays();
|
||||||
|
@ -495,8 +491,8 @@ public class GraggBulirschStoerIntegrator
|
||||||
for (int j = 1; j < k; ++j) {
|
for (int j = 1; j < k; ++j) {
|
||||||
for (int i = 0; i < last.length; ++i) {
|
for (int i = 0; i < last.length; ++i) {
|
||||||
// Aitken-Neville's recursive formula
|
// Aitken-Neville's recursive formula
|
||||||
diag[k-j-1][i] = diag[k-j][i]
|
diag[k-j-1][i] = diag[k-j][i] +
|
||||||
+ coeff[k+offset][j-1] * (diag[k-j][i] - diag[k-j-1][i]);
|
coeff[k+offset][j-1] * (diag[k-j][i] - diag[k-j-1][i]);
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
|
@ -563,11 +559,9 @@ public class GraggBulirschStoerIntegrator
|
||||||
rescale(y, y, scale);
|
rescale(y, y, scale);
|
||||||
|
|
||||||
// initial order selection
|
// initial order selection
|
||||||
double log10R = Math.log(Math.max(1.0e-10,
|
double tol =
|
||||||
(vecRelativeTolerance == null)
|
(vecRelativeTolerance == null) ? scalRelativeTolerance : vecRelativeTolerance[0];
|
||||||
? scalRelativeTolerance
|
double log10R = Math.log(Math.max(1.0e-10, tol)) / Math.log(10.0);
|
||||||
: vecRelativeTolerance[0]))
|
|
||||||
/ Math.log(10.0);
|
|
||||||
int targetIter = Math.max(1,
|
int targetIter = Math.max(1,
|
||||||
Math.min(sequence.length - 2,
|
Math.min(sequence.length - 2,
|
||||||
(int) Math.floor(0.5 - 0.6 * log10R)));
|
(int) Math.floor(0.5 - 0.6 * log10R)));
|
||||||
|
@ -625,8 +619,8 @@ public class GraggBulirschStoerIntegrator
|
||||||
stepSize = hNew;
|
stepSize = hNew;
|
||||||
|
|
||||||
// step adjustment near bounds
|
// step adjustment near bounds
|
||||||
if ((forward && (stepStart + stepSize > t))
|
if ((forward && (stepStart + stepSize > t)) ||
|
||||||
|| ((! forward) && (stepStart + stepSize < t))) {
|
((! forward) && (stepStart + stepSize < t))) {
|
||||||
stepSize = t - stepStart;
|
stepSize = t - stepStart;
|
||||||
}
|
}
|
||||||
double nextT = stepStart + stepSize;
|
double nextT = stepStart + stepSize;
|
||||||
|
@ -698,17 +692,17 @@ public class GraggBulirschStoerIntegrator
|
||||||
// estimate if there is a chance convergence will
|
// estimate if there is a chance convergence will
|
||||||
// be reached on next iteration, using the
|
// be reached on next iteration, using the
|
||||||
// asymptotic evolution of error
|
// asymptotic evolution of error
|
||||||
double ratio = ((double) sequence [k] * sequence[k+1])
|
double ratio = ((double) sequence [k] * sequence[k+1]) /
|
||||||
/ (sequence[0] * sequence[0]);
|
(sequence[0] * sequence[0]);
|
||||||
if (error > ratio * ratio) {
|
if (error > ratio * ratio) {
|
||||||
// we don't expect to converge on next iteration
|
// we don't expect to converge on next iteration
|
||||||
// we reject the step immediately and reduce order
|
// we reject the step immediately and reduce order
|
||||||
reject = true;
|
reject = true;
|
||||||
loop = false;
|
loop = false;
|
||||||
targetIter = k;
|
targetIter = k;
|
||||||
if ((targetIter > 1)
|
if ((targetIter > 1) &&
|
||||||
&& (costPerTimeUnit[targetIter-1]
|
(costPerTimeUnit[targetIter-1] <
|
||||||
< orderControl1 * costPerTimeUnit[targetIter])) {
|
orderControl1 * costPerTimeUnit[targetIter])) {
|
||||||
--targetIter;
|
--targetIter;
|
||||||
}
|
}
|
||||||
hNew = optimalStep[targetIter];
|
hNew = optimalStep[targetIter];
|
||||||
|
@ -731,9 +725,9 @@ public class GraggBulirschStoerIntegrator
|
||||||
// we reject the step immediately
|
// we reject the step immediately
|
||||||
reject = true;
|
reject = true;
|
||||||
loop = false;
|
loop = false;
|
||||||
if ((targetIter > 1)
|
if ((targetIter > 1) &&
|
||||||
&& (costPerTimeUnit[targetIter-1]
|
(costPerTimeUnit[targetIter-1] <
|
||||||
< orderControl1 * costPerTimeUnit[targetIter])) {
|
orderControl1 * costPerTimeUnit[targetIter])) {
|
||||||
--targetIter;
|
--targetIter;
|
||||||
}
|
}
|
||||||
hNew = optimalStep[targetIter];
|
hNew = optimalStep[targetIter];
|
||||||
|
@ -744,9 +738,9 @@ public class GraggBulirschStoerIntegrator
|
||||||
case 1 :
|
case 1 :
|
||||||
if (error > 1.0) {
|
if (error > 1.0) {
|
||||||
reject = true;
|
reject = true;
|
||||||
if ((targetIter > 1)
|
if ((targetIter > 1) &&
|
||||||
&& (costPerTimeUnit[targetIter-1]
|
(costPerTimeUnit[targetIter-1] <
|
||||||
< orderControl1 * costPerTimeUnit[targetIter])) {
|
orderControl1 * costPerTimeUnit[targetIter])) {
|
||||||
--targetIter;
|
--targetIter;
|
||||||
}
|
}
|
||||||
hNew = optimalStep[targetIter];
|
hNew = optimalStep[targetIter];
|
||||||
|
@ -887,8 +881,8 @@ public class GraggBulirschStoerIntegrator
|
||||||
}
|
}
|
||||||
} else {
|
} else {
|
||||||
optimalIter = k - 1;
|
optimalIter = k - 1;
|
||||||
if ((k > 2)
|
if ((k > 2) &&
|
||||||
&& (costPerTimeUnit[k-2] < orderControl1 * costPerTimeUnit[k-1])) {
|
(costPerTimeUnit[k-2] < orderControl1 * costPerTimeUnit[k-1])) {
|
||||||
optimalIter = k - 2;
|
optimalIter = k - 2;
|
||||||
}
|
}
|
||||||
if (costPerTimeUnit[k] < orderControl2 * costPerTimeUnit[optimalIter]) {
|
if (costPerTimeUnit[k] < orderControl2 * costPerTimeUnit[optimalIter]) {
|
||||||
|
@ -906,14 +900,14 @@ public class GraggBulirschStoerIntegrator
|
||||||
if (optimalIter <= k) {
|
if (optimalIter <= k) {
|
||||||
hNew = optimalStep[optimalIter];
|
hNew = optimalStep[optimalIter];
|
||||||
} else {
|
} else {
|
||||||
if ((k < targetIter)
|
if ((k < targetIter) &&
|
||||||
&& (costPerTimeUnit[k] < orderControl2 * costPerTimeUnit[k-1])) {
|
(costPerTimeUnit[k] < orderControl2 * costPerTimeUnit[k-1])) {
|
||||||
hNew = filterStep(optimalStep[k]
|
hNew = filterStep(optimalStep[k] *
|
||||||
* costPerStep[optimalIter+1] / costPerStep[k],
|
costPerStep[optimalIter+1] / costPerStep[k],
|
||||||
false);
|
false);
|
||||||
} else {
|
} else {
|
||||||
hNew = filterStep(optimalStep[k]
|
hNew = filterStep(optimalStep[k] *
|
||||||
* costPerStep[optimalIter] / costPerStep[k],
|
costPerStep[optimalIter] / costPerStep[k],
|
||||||
false);
|
false);
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
|
@ -274,9 +274,8 @@ class GraggBulirschStoerStepInterpolator
|
||||||
for (int j = 4; j <= mu; ++j) {
|
for (int j = 4; j <= mu; ++j) {
|
||||||
double fac1 = 0.5 * j * (j - 1);
|
double fac1 = 0.5 * j * (j - 1);
|
||||||
double fac2 = 2 * fac1 * (j - 2) * (j - 3);
|
double fac2 = 2 * fac1 * (j - 2) * (j - 3);
|
||||||
polynoms[j+4][i] = 16 * (yMidDots[j][i]
|
polynoms[j+4][i] =
|
||||||
+ fac1 * polynoms[j+2][i]
|
16 * (yMidDots[j][i] + fac1 * polynoms[j+2][i] - fac2 * polynoms[j][i]);
|
||||||
- fac2 * polynoms[j][i]);
|
|
||||||
}
|
}
|
||||||
|
|
||||||
}
|
}
|
||||||
|
@ -324,10 +323,10 @@ class GraggBulirschStoerStepInterpolator
|
||||||
t4 = t4 * t4;
|
t4 = t4 * t4;
|
||||||
|
|
||||||
for (int i = 0; i < dimension; ++i) {
|
for (int i = 0; i < dimension; ++i) {
|
||||||
interpolatedState[i] = polynoms[0][i]
|
interpolatedState[i] = polynoms[0][i] +
|
||||||
+ theta * (polynoms[1][i]
|
theta * (polynoms[1][i] +
|
||||||
+ oneMinusTheta * (polynoms[2][i] * theta
|
oneMinusTheta * (polynoms[2][i] * theta +
|
||||||
+ polynoms[3][i] * oneMinusTheta));
|
polynoms[3][i] * oneMinusTheta));
|
||||||
|
|
||||||
if (currentDegree > 3) {
|
if (currentDegree > 3) {
|
||||||
double c = polynoms[currentDegree][i];
|
double c = polynoms[currentDegree][i];
|
||||||
|
|
|
@ -123,9 +123,9 @@ public class HighamHall54Integrator
|
||||||
}
|
}
|
||||||
|
|
||||||
double yScale = Math.max(Math.abs(y0[j]), Math.abs(y1[j]));
|
double yScale = Math.max(Math.abs(y0[j]), Math.abs(y1[j]));
|
||||||
double tol = (vecAbsoluteTolerance == null)
|
double tol = (vecAbsoluteTolerance == null) ?
|
||||||
? (scalAbsoluteTolerance + scalRelativeTolerance * yScale)
|
(scalAbsoluteTolerance + scalRelativeTolerance * yScale) :
|
||||||
: (vecAbsoluteTolerance[j] + vecRelativeTolerance[j] * yScale);
|
(vecAbsoluteTolerance[j] + vecRelativeTolerance[j] * yScale);
|
||||||
double ratio = h * errSum / tol;
|
double ratio = h * errSum / tol;
|
||||||
error += ratio * ratio;
|
error += ratio * ratio;
|
||||||
|
|
||||||
|
|
|
@ -80,9 +80,9 @@ class HighamHall54StepInterpolator
|
||||||
double b5 = h * (-5.0/48.0 + theta2 * (-5.0/16.0 + theta * 5.0/12.0));
|
double b5 = h * (-5.0/48.0 + theta2 * (-5.0/16.0 + theta * 5.0/12.0));
|
||||||
|
|
||||||
for (int i = 0; i < interpolatedState.length; ++i) {
|
for (int i = 0; i < interpolatedState.length; ++i) {
|
||||||
interpolatedState[i] = currentState[i]
|
interpolatedState[i] = currentState[i] +
|
||||||
+ b0 * yDotK[0][i] + b2 * yDotK[2][i] + b3 * yDotK[3][i]
|
b0 * yDotK[0][i] + b2 * yDotK[2][i] + b3 * yDotK[3][i] +
|
||||||
+ b4 * yDotK[4][i] + b5 * yDotK[5][i];
|
b4 * yDotK[4][i] + b5 * yDotK[5][i];
|
||||||
}
|
}
|
||||||
|
|
||||||
}
|
}
|
||||||
|
|
|
@ -88,8 +88,8 @@ class MidpointStepInterpolator
|
||||||
double coeff2 = oneMinusThetaH * (1.0 + theta);
|
double coeff2 = oneMinusThetaH * (1.0 + theta);
|
||||||
|
|
||||||
for (int i = 0; i < interpolatedState.length; ++i) {
|
for (int i = 0; i < interpolatedState.length; ++i) {
|
||||||
interpolatedState[i] = currentState[i]
|
interpolatedState[i] = currentState[i] +
|
||||||
+ coeff1 * yDotK[0][i] - coeff2 * yDotK[1][i];
|
coeff1 * yDotK[0][i] - coeff2 * yDotK[1][i];
|
||||||
}
|
}
|
||||||
|
|
||||||
}
|
}
|
||||||
|
|
|
@ -274,8 +274,8 @@ class SwitchState implements Serializable {
|
||||||
pendingEvent = false;
|
pendingEvent = false;
|
||||||
pendingEventTime = Double.NaN;
|
pendingEventTime = Double.NaN;
|
||||||
|
|
||||||
return (nextAction == SwitchingFunction.RESET_STATE)
|
return (nextAction == SwitchingFunction.RESET_STATE) ||
|
||||||
|| (nextAction == SwitchingFunction.RESET_DERIVATIVES);
|
(nextAction == SwitchingFunction.RESET_DERIVATIVES);
|
||||||
|
|
||||||
}
|
}
|
||||||
|
|
||||||
|
|
|
@ -97,9 +97,9 @@ class ThreeEighthesStepInterpolator
|
||||||
double coeff4 = s * (1 + theta + fourTheta2);
|
double coeff4 = s * (1 + theta + fourTheta2);
|
||||||
|
|
||||||
for (int i = 0; i < interpolatedState.length; ++i) {
|
for (int i = 0; i < interpolatedState.length; ++i) {
|
||||||
interpolatedState[i] = currentState[i]
|
interpolatedState[i] = currentState[i] -
|
||||||
- coeff1 * yDotK[0][i] - coeff2 * yDotK[1][i]
|
coeff1 * yDotK[0][i] - coeff2 * yDotK[1][i] -
|
||||||
- coeff3 * yDotK[2][i] - coeff4 * yDotK[3][i];
|
coeff3 * yDotK[2][i] - coeff4 * yDotK[3][i];
|
||||||
}
|
}
|
||||||
|
|
||||||
}
|
}
|
||||||
|
|
|
@ -504,11 +504,11 @@ public abstract class DirectSearchOptimizer {
|
||||||
|
|
||||||
// return the found point given the lowest cost
|
// return the found point given the lowest cost
|
||||||
if (minima[0] == null) {
|
if (minima[0] == null) {
|
||||||
throw new ConvergenceException("none of the {0} start points"
|
throw new ConvergenceException("none of the {0} start points" +
|
||||||
+ " lead to convergence",
|
" lead to convergence",
|
||||||
new String[] {
|
new Object[] {
|
||||||
Integer.toString(starts)
|
Integer.toString(starts)
|
||||||
});
|
});
|
||||||
}
|
}
|
||||||
return minima[0];
|
return minima[0];
|
||||||
|
|
||||||
|
|
Loading…
Reference in New Issue