List of usage examples for org.apache.commons.math3.ode FirstOrderDifferentialEquations FirstOrderDifferentialEquations
FirstOrderDifferentialEquations
From source file:eu.mihosoft.fx.tutorials.gravity.SolarSystem.java
/** * Creates a system of odes that describe the physical behavior of the solar * system./*from w w w. j a v a 2 s .com*/ * * @param numBodies number of bodies to simulate * @param y state vector (location and velocity) * @param m particle masses * @param ignoreFlags ignore flags (for simulating collisions) * @param G gravitational constant * @return system of first order differential equations */ private FirstOrderDifferentialEquations createODE(final int numBodies, final double[] y, final double[] m, final boolean[] ignoreFlags, final double G) { return new FirstOrderDifferentialEquations() { @Override public int getDimension() { return numBodies * ODEParticle.getStructSize(); } @Override public void computeDerivatives(double t, double[] y, double[] yDot) throws MaxCountExceededException, DimensionMismatchException { ODEParticle pI = new ODEParticle(y, m, ignoreFlags, 0); ODEParticle pJ = new ODEParticle(y, m, ignoreFlags, 1); // d2 rI / dt^2 = GmJ*(rJ-rI)/|rJ-rI|^3 // http://www.physics.buffalo.edu/phy410-505/topic5/ // http://physics.princeton.edu/~fpretori/Nbody/intro.htm // http://de.wikipedia.org/wiki/Newtonsches_Gravitationsgesetz // outer sum for (int i = 0; i < numBodies; i++) { double aIX = 0; double aIY = 0; double aIZ = 0; pI.setIndex(i); if (pI.isIgnored()) { continue; } // inner sum for (int j = 0; j < numBodies; j++) { if (i == j) { continue; } pJ.setIndex(j); if (pJ.isIgnored()) { continue; } double rJMinusRIX = pJ.getRX() - pI.getRX(); double rJMinusRIY = pJ.getRY() - pI.getRY(); double rJMinusRIZ = pJ.getRZ() - pI.getRZ(); double magnitudeRJMinusRISquare = rJMinusRIX * rJMinusRIX + rJMinusRIY * rJMinusRIY + rJMinusRIZ * rJMinusRIZ; double magnitudeRJMinusRI = Math.sqrt(magnitudeRJMinusRISquare); aIX += pI.getMass() * pJ.getMass() * rJMinusRIX / (magnitudeRJMinusRISquare * magnitudeRJMinusRI); aIY += pI.getMass() * pJ.getMass() * rJMinusRIY / (magnitudeRJMinusRISquare * magnitudeRJMinusRI); aIZ += pI.getMass() * pJ.getMass() * rJMinusRIZ / (magnitudeRJMinusRISquare * magnitudeRJMinusRI); } pI.setRDerivativeTo(yDot); aIX *= G / pI.getMass(); aIY *= G / pI.getMass(); aIZ *= G / pI.getMass(); pI.setVDerivativeTo(yDot, aIX, aIY, aIZ); } } }; }
From source file:de.bund.bfr.math.MathUtils.java
public static FirstOrderDifferentialEquations createDiffEquations(Parser parser, List<ASTNode> functions, List<String> dependentVariables, String timeVariable, Map<String, UnivariateFunction> variableFunctions) { return new FirstOrderDifferentialEquations() { @Override// w w w .j a v a2 s. com public int getDimension() { return functions.size(); } @Override public void computeDerivatives(double t, double[] y, double[] yDot) throws MaxCountExceededException, DimensionMismatchException { parser.setVarValue(timeVariable, t); variableFunctions.forEach((var, function) -> parser.setVarValue(var, function.value(t))); for (int i = 0; i < functions.size(); i++) { parser.setVarValue(dependentVariables.get(i), y[i]); } for (int i = 0; i < functions.size(); i++) { try { double value = parser.evaluate(functions.get(i)); yDot[i] = Double.isFinite(value) ? value : Double.NaN; } catch (ParseException e) { e.printStackTrace(); yDot[i] = Double.NaN; } } } }; }
From source file:nl.rivm.cib.episim.model.disease.infection.MSEIRSTest.java
public static Observable<Map.Entry<Double, double[]>> deterministic(final SIRConfig config, final Supplier<FirstOrderIntegrator> integrators) { return Observable.create(sub -> { final double gamma = 1. / config.recovery(); final double beta = gamma * config.reproduction(); final double[] y0 = Arrays.stream(config.population()).mapToDouble(n -> n).toArray(); final double[] t = config.t(); try {/*from ww w . jav a 2 s . com*/ final FirstOrderIntegrator integrator = integrators.get(); integrator.addStepHandler(new StepHandler() { @Override public void init(final double t0, final double[] y0, final double t) { publishCopy(sub, t0, y0); } @Override public void handleStep(final StepInterpolator interpolator, final boolean isLast) throws MaxCountExceededException { publishCopy(sub, interpolator.getInterpolatedTime(), interpolator.getInterpolatedState()); if (isLast) sub.onComplete(); } }); integrator.integrate(new FirstOrderDifferentialEquations() { @Override public int getDimension() { return y0.length; } @Override public void computeDerivatives(final double t, final double[] y, final double[] yp) { // SIR terms (flow rates) final double n = y[0] + y[1] + y[2], flow_si = beta * y[0] * y[1] / n, flow_ir = gamma * y[1]; yp[0] = -flow_si; yp[1] = flow_si - flow_ir; yp[2] = flow_ir; } }, t[0], y0, t[1], y0); } catch (final Exception e) { sub.onError(e); } }); }
From source file:org.orekit.utils.AngularCoordinatesTest.java
@Test public void testShiftWithAcceleration() throws OrekitException { double rate = 2 * FastMath.PI / (12 * 60); double acc = 0.001; double dt = 1.0; int n = 2000; final AngularCoordinates quadratic = new AngularCoordinates(Rotation.IDENTITY, new Vector3D(rate, Vector3D.PLUS_K), new Vector3D(acc, Vector3D.PLUS_J)); final AngularCoordinates linear = new AngularCoordinates(quadratic.getRotation(), quadratic.getRotationRate(), Vector3D.ZERO); final FirstOrderDifferentialEquations ode = new FirstOrderDifferentialEquations() { public int getDimension() { return 4; }// w ww .j av a 2 s . c om public void computeDerivatives(final double t, final double[] q, final double[] qDot) { final double omegaX = quadratic.getRotationRate().getX() + t * quadratic.getRotationAcceleration().getX(); final double omegaY = quadratic.getRotationRate().getY() + t * quadratic.getRotationAcceleration().getY(); final double omegaZ = quadratic.getRotationRate().getZ() + t * quadratic.getRotationAcceleration().getZ(); qDot[0] = 0.5 * MathArrays.linearCombination(-q[1], omegaX, -q[2], omegaY, -q[3], omegaZ); qDot[1] = 0.5 * MathArrays.linearCombination(q[0], omegaX, -q[3], omegaY, q[2], omegaZ); qDot[2] = 0.5 * MathArrays.linearCombination(q[3], omegaX, q[0], omegaY, -q[1], omegaZ); qDot[3] = 0.5 * MathArrays.linearCombination(-q[2], omegaX, q[1], omegaY, q[0], omegaZ); } }; FirstOrderIntegrator integrator = new DormandPrince853Integrator(1.0e-6, 1.0, 1.0e-12, 1.0e-12); integrator.addStepHandler(new StepNormalizer(dt / n, new FixedStepHandler() { public void init(double t0, double[] y0, double t) { } public void handleStep(double t, double[] y, double[] yDot, boolean isLast) { Rotation reference = new Rotation(y[0], y[1], y[2], y[3], true); // the error in shiftedBy taking acceleration into account is cubic double expectedCubicError = 1.4544e-6 * t * t * t; Assert.assertEquals(expectedCubicError, Rotation.distance(reference, quadratic.shiftedBy(t).getRotation()), 0.0001 * expectedCubicError); // the error in shiftedBy not taking acceleration into account is quadratic double expectedQuadraticError = 5.0e-4 * t * t; Assert.assertEquals(expectedQuadraticError, Rotation.distance(reference, linear.shiftedBy(t).getRotation()), 0.00001 * expectedQuadraticError); } })); double[] y = new double[] { quadratic.getRotation().getQ0(), quadratic.getRotation().getQ1(), quadratic.getRotation().getQ2(), quadratic.getRotation().getQ3() }; integrator.integrate(ode, 0, y, dt, y); }
From source file:org.orekit.utils.TimeStampedAngularCoordinatesTest.java
private double[] interpolationErrors(final TimeStampedAngularCoordinates reference, double dt) throws OrekitException { final FirstOrderDifferentialEquations ode = new FirstOrderDifferentialEquations() { public int getDimension() { return 4; }// w w w.j av a2s. co m public void computeDerivatives(final double t, final double[] q, final double[] qDot) { final double omegaX = reference.getRotationRate().getX() + t * reference.getRotationAcceleration().getX(); final double omegaY = reference.getRotationRate().getY() + t * reference.getRotationAcceleration().getY(); final double omegaZ = reference.getRotationRate().getZ() + t * reference.getRotationAcceleration().getZ(); qDot[0] = 0.5 * MathArrays.linearCombination(-q[1], omegaX, -q[2], omegaY, -q[3], omegaZ); qDot[1] = 0.5 * MathArrays.linearCombination(q[0], omegaX, -q[3], omegaY, q[2], omegaZ); qDot[2] = 0.5 * MathArrays.linearCombination(q[3], omegaX, q[0], omegaY, -q[1], omegaZ); qDot[3] = 0.5 * MathArrays.linearCombination(-q[2], omegaX, q[1], omegaY, q[0], omegaZ); } }; final List<TimeStampedAngularCoordinates> complete = new ArrayList<TimeStampedAngularCoordinates>(); FirstOrderIntegrator integrator = new DormandPrince853Integrator(1.0e-6, 1.0, 1.0e-12, 1.0e-12); integrator.addStepHandler(new StepNormalizer(dt / 2000, new FixedStepHandler() { public void init(double t0, double[] y0, double t) { } public void handleStep(double t, double[] y, double[] yDot, boolean isLast) { complete.add( new TimeStampedAngularCoordinates(reference.getDate().shiftedBy(t), new Rotation(y[0], y[1], y[2], y[3], true), new Vector3D(1, reference.getRotationRate(), t, reference.getRotationAcceleration()), reference.getRotationAcceleration())); } })); double[] y = new double[] { reference.getRotation().getQ0(), reference.getRotation().getQ1(), reference.getRotation().getQ2(), reference.getRotation().getQ3() }; integrator.integrate(ode, 0, y, dt, y); List<TimeStampedAngularCoordinates> sample = new ArrayList<TimeStampedAngularCoordinates>(); sample.add(complete.get(0)); sample.add(complete.get(complete.size() / 2)); sample.add(complete.get(complete.size() - 1)); double maxRotationError = 0; double maxRateError = 0; double maxAccelerationError = 0; for (TimeStampedAngularCoordinates acRef : complete) { TimeStampedAngularCoordinates interpolated = TimeStampedAngularCoordinates.interpolate(acRef.getDate(), AngularDerivativesFilter.USE_RRA, sample); double rotationError = Rotation.distance(acRef.getRotation(), interpolated.getRotation()); double rateError = Vector3D.distance(acRef.getRotationRate(), interpolated.getRotationRate()); double accelerationError = Vector3D.distance(acRef.getRotationAcceleration(), interpolated.getRotationAcceleration()); maxRotationError = FastMath.max(maxRotationError, rotationError); maxRateError = FastMath.max(maxRateError, rateError); maxAccelerationError = FastMath.max(maxAccelerationError, accelerationError); } return new double[] { maxRotationError, maxRateError, maxAccelerationError }; }