List of usage examples for org.apache.commons.math3.ode FirstOrderIntegrator integrate
double integrate(FirstOrderDifferentialEquations equations, double t0, double[] y0, double t, double[] y) throws DimensionMismatchException, NumberIsTooSmallException, MaxCountExceededException, NoBracketingException;
From source file:jat.examples.TwoBodyExample.TwoBodyExample.java
public static void main(String[] args) { TwoBodyExample x = new TwoBodyExample(); // create a TwoBody orbit using orbit elements TwoBodyAPL sat = new TwoBodyAPL(7000.0, 0.3, 0.0, 0.0, 0.0, 0.0); double[] y = sat.randv(); ArrayRealVector v = new ArrayRealVector(y); DecimalFormat df2 = new DecimalFormat("#,###,###,##0.00"); RealVectorFormat format = new RealVectorFormat(df2); System.out.println(format.format(v)); // find out the period of the orbit double period = sat.period(); // set the final time = one orbit period double tf = period; // set the initial time to zero double t0 = 0.0; // propagate the orbit FirstOrderIntegrator dp853 = new DormandPrince853Integrator(1.0e-8, 100.0, 1.0e-10, 1.0e-10); dp853.addStepHandler(sat.stepHandler); // double[] y = new double[] { 7000.0, 0, 0, .0, 8, 0 }; // initial // state/*w w w . j a v a 2 s . c o m*/ dp853.integrate(sat, 0.0, y, 8000, y); // now y contains final state at // tf Double[] objArray = sat.time.toArray(new Double[sat.time.size()]); double[] timeArray = ArrayUtils.toPrimitive(objArray); double[] xsolArray = ArrayUtils.toPrimitive(sat.xsol.toArray(new Double[sat.time.size()])); double[] ysolArray = ArrayUtils.toPrimitive(sat.ysol.toArray(new Double[sat.time.size()])); double[][] XY = new double[timeArray.length][2]; // int a=0; // System.arraycopy(timeArray,0,XY[a],0,timeArray.length); // System.arraycopy(ysolArray,0,XY[1],0,ysolArray.length); for (int i = 0; i < timeArray.length; i++) { XY[i][0] = xsolArray[i]; XY[i][1] = ysolArray[i]; } Plot2DPanel p = new Plot2DPanel(); // Plot2DPanel p = new Plot2DPanel(min, max, axesScales, axesLabels); ScatterPlot s = new ScatterPlot("orbit", Color.RED, XY); // LinePlot l = new LinePlot("sin", Color.RED, XY); // l.closed_curve = false; // l.draw_dot = true; p.addPlot(s); p.setLegendOrientation(PlotPanel.SOUTH); double plotSize = 10000.; double[] min = { -plotSize, -plotSize }; double[] max = { plotSize, plotSize }; p.setFixedBounds(min, max); new FrameView(p).setDefaultCloseOperation(JFrame.EXIT_ON_CLOSE); System.out.println("end"); }
From source file:de.bund.bfr.math.Evaluator.java
public static double[] getDiffPoints(Map<String, Double> parserConstants, Map<String, String> functions, Map<String, Double> initValues, Map<String, String> initParameters, Map<String, List<Double>> conditionLists, String dependentVariable, Map<String, Double> independentVariables, String varX, double[] valuesX, IntegratorFactory integrator, InterpolationFactory interpolator) throws ParseException { DiffFunctionConf function = new DiffFunctionConf(parserConstants, functions, initValues, initParameters, conditionLists, dependentVariable, independentVariables, varX, valuesX, integrator, interpolator); double[] result = diffResults.getIfPresent(function); if (result != null) { return result; }/* w w w .j ava2s. co m*/ List<ASTNode> fs = new ArrayList<>(); List<String> valueVariables = new ArrayList<>(); double[] values = new double[functions.size()]; Parser parser = new Parser(); parserConstants.forEach((constant, value) -> parser.setVarValue(constant, value)); int index = 0; for (Map.Entry<String, String> entry : functions.entrySet()) { String var = entry.getKey(); fs.add(parser.parse(entry.getValue())); valueVariables.add(var); values[index++] = initValues.containsKey(var) ? initValues.get(var) : parserConstants.get(initParameters.get(var)); } Map<String, UnivariateFunction> variableFunctions = MathUtils.createInterpolationFunctions(conditionLists, varX, interpolator); FirstOrderDifferentialEquations f = MathUtils.createDiffEquations(parser, fs, valueVariables, varX, variableFunctions); FirstOrderIntegrator instance = integrator.createIntegrator(); double diffValue = conditionLists.get(varX).get(0); int depIndex = valueVariables.indexOf(dependentVariable); double[] valuesY = new double[valuesX.length]; for (int i = 0; i < valuesX.length; i++) { if (valuesX[i] == diffValue) { valuesY[i] = values[depIndex]; } else if (valuesX[i] > diffValue) { instance.integrate(f, diffValue, values, valuesX[i], values); diffValue = valuesX[i]; valuesY[i] = values[depIndex]; } else { valuesY[i] = Double.NaN; } } diffResults.put(function, valuesY); return valuesY; }
From source file:jat.examples.ephemeris.DE405PropagateText.java
void doExample() { //double tf = 10000000.; double tf = 3600 * 24 * 500; double[] y = { 2e8, 0, 0, 0, 24.2, 0 }; // initial state PathUtil path = new PathUtil(); DE405Plus Eph = new DE405Plus(path); Eph.setFrame(DE405Frame.frame.HEE); Eph.printSteps = true;//from w ww . ja va 2 s.c o m TimeAPL myTime = new TimeAPL(2003, 3, 1, 12, 0, 0); Eph.setIntegrationStartTime(myTime); Eph.bodyGravOnOff[body.SUN.ordinal()] = true; FirstOrderIntegrator dp853 = new DormandPrince853Integrator(1.0e-8, tf / 10.0, 1.0e-10, 1.0e-10); dp853.addStepHandler(Eph.stepHandler); FirstOrderDifferentialEquations ode = Eph; dp853.integrate(ode, 0.0, y, tf, y); // now y contains final state // at // time tf if (print) { String nf = "%10.3f "; String format = nf + nf + nf + nf + nf; System.out.printf(format, tf, y[0], y[1], y[2], Eph.energy(tf, y)); System.out.println(); } }
From source file:jat.examples.ephemeris.DE405PropagatePlot.java
void doExample() { double tf = 3600 * 24 * 300; double[] y0 = { 2e8, 0, 0, 0, 24.2, 0 }; // initial state double[] y = new double[6]; PathUtil path = new PathUtil(); DE405Plus Eph = new DE405Plus(path); Eph.setFrame(DE405Frame.frame.HEE); Eph.printSteps = true;// ww w . ja v a 2 s. co m TimeAPL myTime = new TimeAPL(2003, 3, 1, 12, 0, 0); Eph.setIntegrationStartTime(myTime); Eph.bodyGravOnOff[body.SUN.ordinal()] = true; // Eph.planetOnOff[body.JUPITER.ordinal()] = true; FirstOrderIntegrator dp853 = new DormandPrince853Integrator(1.0e-8, tf / 10.0, 1.0e-10, 1.0e-10); dp853.addStepHandler(Eph.stepHandler); FirstOrderDifferentialEquations ode = Eph; dp853.integrate(ode, 0.0, y0, tf, y); // now y contains final state at // time tf if (print) { String nf = "%10.3f "; String format = nf + nf + nf + nf + nf; System.out.printf(format, tf, y[0], y[1], y[2], Eph.energy(tf, y)); System.out.println(); } Plot2DPanel p = new Plot2DPanel(); LinePlot l1 = new LinePlot("Jup. off", Color.RED, getXYforPlot(Eph.xsol, Eph.ysol)); l1.closed_curve = false; p.addPlot(l1); Eph.reset(); Eph.bodyGravOnOff[body.JUPITER.ordinal()] = true; dp853.integrate(ode, 0.0, y0, tf, y); // now y contains final state at LinePlot l2 = new LinePlot("Jup. on", Color.BLUE, getXYforPlot(Eph.xsol, Eph.ysol)); l2.closed_curve = false; p.addPlot(l2); VectorN EarthPos = null; try { EarthPos = Eph.get_planet_pos(body.EARTH, myTime); } catch (IOException e) { // TODO Auto-generated catch block e.printStackTrace(); } addPoint(p, "Earth", java.awt.Color.BLUE, EarthPos.x[0], EarthPos.x[1]); p.setLegendOrientation(PlotPanel.SOUTH); double plotSize = 2e8; p.setFixedBounds(0, -plotSize, plotSize); p.setFixedBounds(1, -plotSize, plotSize); new FrameView(p).setDefaultCloseOperation(JFrame.EXIT_ON_CLOSE); }
From source file:de.bund.bfr.math.VectorDiffFunction.java
@Override public double[] value(double[] point) throws IllegalArgumentException { for (int i = 0; i < parameters.size(); i++) { if (!initParameters.contains(parameters.get(i))) { parser.setVarValue(parameters.get(i), point[i]); }/*from w w w. j a va2 s. c om*/ } double[] values = new double[formulas.size()]; for (int i = 0; i < formulas.size(); i++) { values[i] = initValues.get(i) != null ? initValues.get(i) : point[parameters.indexOf(initParameters.get(i))]; } FirstOrderDifferentialEquations f = MathUtils.createDiffEquations(parser, functions, dependentVariables, timeVariable, variableFunctions); FirstOrderIntegrator integratorInstance = integrator.createIntegrator(); List<Double> result = new ArrayList<>(); result.add(values[dependentIndex]); for (int i = 1; i < timeValues.size(); i++) { integratorInstance.integrate(f, timeValues.get(i - 1), values, timeValues.get(i), values); result.add(values[dependentIndex]); } return Doubles.toArray(result); }
From source file:jat.application.DE405Propagator.DE405PropagatorPlot.java
public void add_scene() { // Update Ephemeris to current user parameters for (body b : body.values()) { Eph.bodyGravOnOff[b.ordinal()] = dpMain.dpParam.bodyGravOnOff[b.ordinal()]; }//www .j a v a 2 s .c om Eph.setIntegrationStartTime(dpMain.dpParam.simulationDate); try { Eph.setEarthMoonPlaneNormal(dpMain.dpParam.simulationDate); } catch (IOException e1) { e1.printStackTrace(); } Eph.setFrame(dpMain.dpParam.Frame); Eph.reset(); // Spacecraft Trajectory FirstOrderIntegrator dp853 = new DormandPrince853Integrator(1.0e-8, dpMain.dpParam.tf, 1.0e-10, 1.0e-10); dp853.addStepHandler(Eph.stepHandler); FirstOrderDifferentialEquations ode = Eph; double[] y = new double[6]; dp853.integrate(ode, 0.0, dpMain.dpParam.y0, dpMain.dpParam.tf, y); if (print) { String nf = "%10.3f "; String format = nf + nf + nf + nf + nf; System.out.printf(format, dpMain.dpParam.tf, y[0], y[1], y[2], Eph.energy(dpMain.dpParam.tf, y)); System.out.println(); } LinePlot l1 = new LinePlot("spacecraft", Color.RED, getXYZforPlot(Eph.xsol, Eph.ysol, Eph.zsol)); l1.closed_curve = false; plot.addPlot(l1); addBodies(); // Vector3D y0v=new Vector3D(dpParam.y0[0],dpParam.y0[1],dpParam.y0[2]); // double plotBounds = 2*y0v.getNorm(); plot.setFixedBounds(0, -plotBounds, plotBounds); plot.setFixedBounds(1, -plotBounds, plotBounds); plot.setFixedBounds(2, -plotBounds, plotBounds); plot.setLegendOrientation(PlotPanel.SOUTH); }
From source file:jat.examples.CRTBP.CRTBPPlot.java
void doExample() { double mu = 0.15; double[] y0 = { .11, 0, 0, 1.35, 1.33, 0 }; // initial state // double mu = 0.1; // double mu = 3.035909999e-6; // double mu = 0.012277471; // double[] y0 = { .1, 0, 0, 2.69, 2.69, 0 }; // initial state // double mu = 0.2; CRTBP myCRTBP = new CRTBP(mu); FirstOrderIntegrator dp853 = new DormandPrince853Integrator(1.0e-8, 100.0, 1.0e-10, 1.0e-10); dp853.addStepHandler(myCRTBP.stepHandler); FirstOrderDifferentialEquations ode = myCRTBP; double tf;/*from www. j a va 2 s.co m*/ double[] y = new double[6]; // initial state // for (int i = 1; i < 2; i++) { // tf = i * 20.; tf = 40.; System.arraycopy(y0, 0, y, 0, 6); dp853.integrate(ode, 0.0, y, tf, y); // now y contains final state // at // time tf if (print) { System.out.printf("%9.6f %9.6f %9.6f %9.6f %9.6f", tf, y[0], y[1], y[2], myCRTBP.JacobiIntegral(y)); System.out.println(); } int arraySize = myCRTBP.time.size(); double[] timeArray = ArrayUtils.toPrimitive(myCRTBP.time.toArray(new Double[arraySize])); double[] xsolArray = ArrayUtils.toPrimitive(myCRTBP.xsol.toArray(new Double[arraySize])); double[] ysolArray = ArrayUtils.toPrimitive(myCRTBP.ysol.toArray(new Double[arraySize])); double[][] XY = new double[timeArray.length][2]; for (int i = 0; i < timeArray.length; i++) { XY[i][0] = xsolArray[i]; XY[i][1] = ysolArray[i]; } Plot2DPanel p = new Plot2DPanel(); LinePlot l = new LinePlot("spacecraft", Color.RED, XY); l.closed_curve = false; l.draw_dot = true; p.addPlot(l); double plotSize = 1.2; myCRTBP.findLibrationPoints(); Color darkGreen = new java.awt.Color(0, 190, 0); addPoint(p, "Earth", java.awt.Color.BLUE, -mu, 0); addPoint(p, "Moon", java.awt.Color.gray, 1 - mu, 0); addPoint(p, "L1", darkGreen, myCRTBP.LibPoints[0].getX(), 0); addPoint(p, "L2", darkGreen, myCRTBP.LibPoints[1].getX(), 0); addPoint(p, "L3", darkGreen, myCRTBP.LibPoints[2].getX(), 0); String Labelmu = "mu = " + myCRTBP.mu; p.addLabel(Labelmu, java.awt.Color.black, 1, .9 * plotSize); String initial = "initial x,v = (" + y0[0] + "," + y0[1] + "),(" + y0[3] + "," + y0[4] + ")"; p.addLabel(initial, java.awt.Color.black, 1, .8 * plotSize); String Jacobi = "spacecraft C = " + myCRTBP.C; p.addLabel(Jacobi, java.awt.Color.black, 1, .7 * plotSize); String L1C = "L1 C = " + myCRTBP.C1; p.addLabel(L1C, java.awt.Color.black, 1, .6 * plotSize); myCRTBP.findZeroVelocity(); int size = myCRTBP.xzv.size(); double[] xzvArray = ArrayUtils.toPrimitive(myCRTBP.xzv.toArray(new Double[size])); double[] yzvArray = ArrayUtils.toPrimitive(myCRTBP.yzv.toArray(new Double[size])); double[][] XYzv = new double[size][2]; for (int i = 0; i < size; i++) { XYzv[i][0] = xzvArray[i]; XYzv[i][1] = yzvArray[i]; } LinePlot lzv = new LinePlot("zero vel", Color.blue, XYzv); lzv.closed_curve = false; lzv.draw_dot = true; p.addPlot(lzv); p.setLegendOrientation(PlotPanel.SOUTH); p.setFixedBounds(0, -plotSize, plotSize); p.setFixedBounds(1, -plotSize, plotSize); new FrameView(p).setDefaultCloseOperation(JFrame.EXIT_ON_CLOSE); }
From source file:edu.gcsc.vrl.commons.math.ode.ODESolver.java
public Trajectory solve(@ParamInfo(name = "Label", options = "value=\"Label 1\"") String label, @ParamInfo(name = "x0", options = "value=0.0D") double t0, @ParamInfo(name = "xn", options = "value=3.0D") double tn, @ParamInfo(name = "y0", options = "value=0.0") double y0, @ParamInfo(name = "Min Step", options = "value=1e-6D") double minStep, @ParamInfo(name = "Max Step", options = "value=1e-2D") double maxStep, @ParamInfo(name = "Abs.Tol.", options = "value=1e-10") double absTol, @ParamInfo(name = "Rel.Tol.", options = "value=1e-10") double relTol, @ParamInfo(name = "RHS") FirstOrderDifferentialEquations rhs) { FirstOrderIntegrator integrator = new DormandPrince853Integrator(minStep, maxStep, absTol, relTol); final Trajectory result = new Trajectory(label); StepHandler stepHandler = new StepHandler() { @Override//w w w. ja v a 2 s .c om public void init(double t0, double[] y0, double t) { result.add(t, y0[0]); } @Override public void handleStep(StepInterpolator interpolator, boolean isLast) { double t = interpolator.getCurrentTime(); double[] y = interpolator.getInterpolatedState(); result.add(t, y[0]); } }; integrator.addStepHandler(stepHandler); double[] y = new double[] { y0 }; // initial state integrator.integrate(rhs, t0, y, tn, y); return result; }
From source file:de.bund.bfr.math.MultiVectorDiffFunction.java
@Override public double[] value(double[] point) throws IllegalArgumentException { for (int i = 0; i < parameters.size(); i++) { if (!initParameters.contains(parameters.get(i))) { parser.setVarValue(parameters.get(i), point[i]); }//from ww w . j a v a 2s . co m } FirstOrderIntegrator integratorInstance = integrator.createIntegrator(); List<Double> result = new ArrayList<>(); for (int i = 0; i < timeValues.size(); i++) { FirstOrderDifferentialEquations f = MathUtils.createDiffEquations(parser, functions, dependentVariables, timeVariable, variableFunctions.get(i)); double[] values = new double[formulas.size()]; for (int j = 0; j < formulas.size(); j++) { values[j] = initValues.get(j) != null ? initValues.get(j) : point[parameters.indexOf(initParameters.get(i).get(j))]; } result.add(values[dependentIndex]); for (int j = 1; j < timeValues.get(i).size(); j++) { integratorInstance.integrate(f, timeValues.get(i).get(j - 1), values, timeValues.get(i).get(j), values); result.add(values[dependentIndex]); } } return Doubles.toArray(result); }
From source file:beast.structuredCoalescent.distribution.IndependentStructuredCoalescent.java
public double calculateLogP() { // newly calculate tree intervals treeIntervalsInput.get().calculateIntervals(); // correctly calculate the daughter nodes at coalescent intervals in the case of // bifurcation or in case two nodes are at the same height treeIntervalsInput.get().swap();//from w w w. j av a 2 s . c o m // Set up ArrayLists for the indices of active lineages and the lineage state probabilities activeLineages = new ArrayList<Integer>(); lineStateProbs = new ArrayList<Double>(); // Compute likelihood at each integration time and tree event starting at final sampling time and moving backwards logP = 0; // set the current time double currTime = 0.0; // total number of intervals final int intervalCount = treeIntervalsInput.get().getIntervalCount(); // interval time counter int t = 0; // initialize the number of lineages nr_lineages = 0; // Captures the probabilities of lineages being in a state double[] p; // Initialize the migration rates matrix double[][] migration_rates = new double[states][states]; int c = 0; for (int k = 0; k < states; k++) { for (int l = 0; l < states; l++) { if (k != l) { migration_rates[k][l] = migrationRatesInput.get().getArrayValue(c); c++; } else { // diagonal migration_rates[k][l] = 0.0; } } } // Initialize the coalescent rates double[] coalescent_rates = new double[states]; for (int k = 0; k < states; k++) { coalescent_rates[k] = coalescentRatesInput.get().getArrayValue(k) / 2;//(epiModelInput.get().getF(t,k,k) / (Y.get(k)*Y.get(k))); } // integrate until there are no more tree intervals do { double nextIntervalTime = treeIntervalsInput.get().getInterval(t); // Length of the current interval final double duration = nextIntervalTime;// - currTime; // if the current interval has a length greater than 0, integrate if (duration > 0) { if (dependentHistory) p = new double[lineStateProbs.size()]; // Captures the probabilities of lineages being in a state else p = new double[lineStateProbs.size() + 1]; // Captures the probabilities of lineages being in a state, last one keeps track of the probability // convert the array list to double[] for (int i = 0; i < lineStateProbs.size(); i++) p[i] = lineStateProbs.get(i); // not needed if (!dependentHistory) p[lineStateProbs.size()] = 1; double[] p_for_ode = new double[p.length]; double ts = 0.0; // If proportial time step is true, set the integration time for the given interval // inverse proportional to the number of lineages if (propTimeStep) ts = timeStep / lineStateProbs.size(); else ts = timeStep; // Never choose a longer time step than the integration window if (duration < (ts / 2)) ts = duration / 2; FirstOrderIntegrator integrator = new ClassicalRungeKuttaIntegrator(ts); // set the odes FirstOrderDifferentialEquations ode = new independent_ode_integrator(migration_rates, coalescent_rates, nr_lineages, states); // integrate integrator.integrate(ode, 0, p, duration, p_for_ode); // If the Dimension is larger than the maximum integer, at least one state prob is below 0 and the step is rejected if (ode.getDimension() == Integer.MAX_VALUE) { System.out.println(lineStateProbs.size()); System.out.println("lalalallal"); return Double.NEGATIVE_INFINITY; } for (int i = 0; i < lineStateProbs.size(); i++) lineStateProbs.set(i, p_for_ode[i]); } // update the time currTime = nextIntervalTime; // event is coalescent event if (treeIntervalsInput.get().getIntervalType(t) == IntervalType.COALESCENT) { logP += coalesce(t); nr_lineages--; } // event is sampling event if (treeIntervalsInput.get().getIntervalType(t) == IntervalType.SAMPLE) { addLineages(t); nr_lineages++; } // update the interval number t++; } while (t < intervalCount); //Compute likelihood of remaining tree intervals (coal events occuring before origin) if (Double.isInfinite(logP)) logP = Double.NEGATIVE_INFINITY; if (max_posterior < logP && logP < 0) { max_posterior = logP; max_mig = new double[states * (states - 1)]; max_coal = new double[states]; for (int i = 0; i < 1; i++) max_mig[i] = migrationRatesInput.get().getArrayValue(i); for (int i = 0; i < 1; i++) max_coal[i] = coalescentRatesInput.get().getArrayValue(i); } return logP; }