diff --git a/src/test/java/org/apache/commons/math3/ode/TestFieldProblem1.java b/src/test/java/org/apache/commons/math3/ode/TestFieldProblem1.java new file mode 100644 index 000000000..f79fad028 --- /dev/null +++ b/src/test/java/org/apache/commons/math3/ode/TestFieldProblem1.java @@ -0,0 +1,78 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +package org.apache.commons.math3.ode; + +import org.apache.commons.math3.Field; +import org.apache.commons.math3.RealFieldElement; +import org.apache.commons.math3.util.MathArrays; + +/** + * This class is used in the junit tests for the ODE integrators. + + *

This specific problem is the following differential equation : + * 

+ *    y' = -y
+ *
+ * the solution of this equation is a simple exponential function : + * 
+ *   y (t) = y (t0) exp (t0-t)
+ *
+ *

+ + * @param the type of the field elements + */ +public class TestFieldProblem1> + extends TestFieldProblemAbstract { + + /** + * Simple constructor. + * @param field field to which elements belong + */ + public TestFieldProblem1(Field field) { + super(field); + setInitialConditions(convert(0.0), convert(1.0, 0.1)); + setFinalConditions(convert(4.0)); + setErrorScale(convert(1.0, 1.0)); + } + + @Override + public T[] doComputeDerivatives(T t, T[] y) { + + final T[] yDot = MathArrays.buildArray(getField(), getDimension()); + + // compute the derivatives + for (int i = 0; i < getDimension(); ++i) { + yDot[i] = y[i].negate(); + } + + return yDot; + + } + + @Override + public T[] computeTheoreticalState(T t) { + final T[] y0 = getInitialState(); + final T[] y = MathArrays.buildArray(getField(), getDimension()); + T c = getInitialTime().subtract(t).exp(); + for (int i = 0; i < getDimension(); ++i) { + y[i] = c.multiply(y0[i]); + } + return y; + } + +} diff --git a/src/test/java/org/apache/commons/math3/ode/TestFieldProblem2.java b/src/test/java/org/apache/commons/math3/ode/TestFieldProblem2.java new file mode 100644 index 000000000..7f5b0e375 --- /dev/null +++ b/src/test/java/org/apache/commons/math3/ode/TestFieldProblem2.java @@ -0,0 +1,78 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +package org.apache.commons.math3.ode; + +import org.apache.commons.math3.Field; +import org.apache.commons.math3.RealFieldElement; +import org.apache.commons.math3.util.MathArrays; + +/** + * This class is used in the junit tests for the ODE integrators. + + *

This specific problem is the following differential equation : + * 

+ *    y' = t^3 - t y
+ *
+ * with the initial condition y (0) = 0. The solution of this equation + * is the following function : + * 
+ *   y (t) = t^2 + 2 (exp (- t^2 / 2) - 1)
+ *
+ *

+ + * @param the type of the field elements + */ +public class TestFieldProblem2> + extends TestFieldProblemAbstract { + + /** + * Simple constructor. + * @param field field to which elements belong + */ + public TestFieldProblem2(Field field) { + super(field); + setInitialConditions(convert(0.0), convert(new double[] { 0.0 })); + setFinalConditions(convert(1.0)); + setErrorScale(convert(new double[] { 1.0 })); + } + + @Override + public T[] doComputeDerivatives(T t, T[] y) { + + final T[] yDot = MathArrays.buildArray(getField(), getDimension()); + // compute the derivatives + for (int i = 0; i < getDimension(); ++i) { + yDot[i] = t.multiply(t.multiply(t).subtract(y[i])); + } + + return yDot; + + } + + @Override + public T[] computeTheoreticalState(T t) { + final T[] y = MathArrays.buildArray(getField(), getDimension()); + T t2 = t.multiply(t); + T c = t2.add(t2.multiply(-0.5).exp().subtract(1).multiply(2)); + for (int i = 0; i < getDimension(); ++i) { + y[i] = c; + } + return y; + } + +} diff --git a/src/test/java/org/apache/commons/math3/ode/TestFieldProblem3.java b/src/test/java/org/apache/commons/math3/ode/TestFieldProblem3.java new file mode 100644 index 000000000..9540e81bc --- /dev/null +++ b/src/test/java/org/apache/commons/math3/ode/TestFieldProblem3.java @@ -0,0 +1,124 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +package org.apache.commons.math3.ode; + +import org.apache.commons.math3.Field; +import org.apache.commons.math3.RealFieldElement; +import org.apache.commons.math3.util.MathArrays; + +/** + * This class is used in the junit tests for the ODE integrators. + + *

This specific problem is the following differential equation : + * 

+ *    y1'' = -y1/r^3  y1 (0) = 1-e  y1' (0) = 0
+ *    y2'' = -y2/r^3  y2 (0) = 0    y2' (0) =sqrt((1+e)/(1-e))
+ *    r = sqrt (y1^2 + y2^2), e = 0.9
+ *
+ * This is a two-body problem in the plane which can be solved by + * Kepler's equation + * 
+ *   y1 (t) = ...
+ *
+ *

+ + * @param the type of the field elements + */ +public class TestFieldProblem3> +extends TestFieldProblemAbstract { + + /** Eccentricity */ + T e; + + /** + * Simple constructor. + * @param field field to which elements belong + * @param e eccentricity + */ + public TestFieldProblem3(Field field, T e) { + super(field); + this.e = e; + T[] y0 = MathArrays.buildArray(field, 4); + y0[0] = e.subtract(1).negate(); + y0[1] = field.getZero(); + y0[2] = field.getZero(); + y0[3] = e.add(1).divide(y0[0]).sqrt(); + setInitialConditions(convert(0.0), y0); + setFinalConditions(convert(20.0)); + setErrorScale(convert(1.0, 1.0, 1.0, 1.0)); + } + + /** + * Simple constructor. + * @param field field to which elements belong + */ + public TestFieldProblem3(Field field) { + this(field, field.getZero().add(0.1)); + } + + @Override + public T[] doComputeDerivatives(T t, T[] y) { + + final T[] yDot = MathArrays.buildArray(getField(), getDimension()); + + // current radius + T r2 = y[0].multiply(y[0]).add(y[1].multiply(y[1])); + T invR3 = r2.multiply(r2.sqrt()).reciprocal(); + + // compute the derivatives + yDot[0] = y[2]; + yDot[1] = y[3]; + yDot[2] = invR3.negate().multiply(y[0]); + yDot[3] = invR3.negate().multiply(y[1]); + + return yDot; + + } + + @Override + public T[] computeTheoreticalState(T t) { + + final T[] y = MathArrays.buildArray(getField(), getDimension()); + + // solve Kepler's equation + T E = t; + T d = convert(0); + T corr = convert(999.0); + for (int i = 0; (i < 50) && (corr.abs().getReal() > 1.0e-12); ++i) { + T f2 = e.multiply(E.sin()); + T f0 = d.subtract(f2); + T f1 = e.multiply(E.cos()).subtract(1).negate(); + T f12 = f1.add(f1); + corr = f0.multiply(f12).divide(f1.multiply(f12).subtract(f0.multiply(f2))); + d = d.subtract(corr); + E = t.add(d); + } + + T cosE = E.cos(); + T sinE = E.sin(); + + y[0] = cosE.subtract(e); + y[1] = e.multiply(e).subtract(1).negate().sqrt().multiply(sinE); + y[2] = sinE.divide(e.multiply(cosE).subtract(1)); + y[3] = e.multiply(e).subtract(1).negate().sqrt().multiply(cosE).divide(e.multiply(cosE).subtract(1).negate()); + + return y; + + } + +} diff --git a/src/test/java/org/apache/commons/math3/ode/TestFieldProblem4.java b/src/test/java/org/apache/commons/math3/ode/TestFieldProblem4.java new file mode 100644 index 000000000..0f76c2506 --- /dev/null +++ b/src/test/java/org/apache/commons/math3/ode/TestFieldProblem4.java @@ -0,0 +1,160 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +package org.apache.commons.math3.ode; + +import java.lang.reflect.Array; + +import org.apache.commons.math3.Field; +import org.apache.commons.math3.RealFieldElement; +import org.apache.commons.math3.ode.events.Action; +import org.apache.commons.math3.ode.events.FieldEventHandler; +import org.apache.commons.math3.util.FastMath; +import org.apache.commons.math3.util.MathArrays; + +/** + * This class is used in the junit tests for the ODE integrators. + + *

This specific problem is the following differential equation : + * 

+ *    x'' = -x
+ *
+ * And when x decreases down to 0, the state should be changed as follows : + * 
+ *   x' -> -x'
+ *
+ * The theoretical solution of this problem is x = |sin(t+a)| + *

+ + * @param the type of the field elements + */ +public class TestFieldProblem4> + extends TestFieldProblemAbstract { + + /** Time offset. */ + private T a; + + /** Simple constructor. + * @param field field to which elements belong + */ + public TestFieldProblem4(Field field) { + super(field); + a = convert(1.2); + T[] y0 = MathArrays.buildArray(field, 2); + y0[0] = a.sin(); + y0[1] = a.cos();; + setInitialConditions(convert(0.0), y0); + setFinalConditions(convert(15)); + setErrorScale(convert(1.0, 0.0)); + } + + @Override + public FieldEventHandler[] getEventsHandlers() { + @SuppressWarnings("unchecked") + FieldEventHandler[] handlers = + (FieldEventHandler[]) Array.newInstance(FieldEventHandler.class, 2); + handlers[0] = new Bounce(); + handlers[1] = new Stop(); + return handlers; + } + + /** + * Get the theoretical events times. + * @return theoretical events times + */ + @Override + public T[] getTheoreticalEventsTimes() { + T[] array = MathArrays.buildArray(getField(), 5); + array[0] = a.negate().add(1 * FastMath.PI); + array[1] = a.negate().add(2 * FastMath.PI); + array[2] = a.negate().add(3 * FastMath.PI); + array[3] = a.negate().add(4 * FastMath.PI); + array[4] = convert(120.0); + return array; + } + + @Override + public T[] doComputeDerivatives(T t, T[] y) { + final T[] yDot = MathArrays.buildArray(getField(), getDimension()); + yDot[0] = y[1]; + yDot[1] = y[0].negate(); + return yDot; + } + + @Override + public T[] computeTheoreticalState(T t) { + T sin = t.add(a).sin(); + T cos = t.add(a).cos(); + final T[] y = MathArrays.buildArray(getField(), getDimension()); + y[0] = sin.abs(); + y[1] = (sin.getReal() >= 0) ? cos : cos.negate(); + return y; + } + + private static class Bounce> implements FieldEventHandler { + + private int sign; + + public Bounce() { + sign = +1; + } + + public void init(FieldODEStateAndDerivative state0, T t) { + } + + public T g(FieldODEStateAndDerivative state) { + return state.getState()[0].multiply(sign); + } + + public Action eventOccurred(FieldODEStateAndDerivative state, boolean increasing) { + // this sign change is needed because the state will be reset soon + sign = -sign; + return Action.RESET_STATE; + } + + public FieldODEState resetState(FieldODEStateAndDerivative state) { + T[] y = state.getState(); + y[0] = y[0].negate(); + y[1] = y[1].negate(); + return new FieldODEState(state.getTime(), y); + } + + } + + private static class Stop> implements FieldEventHandler { + + public Stop() { + } + + public void init(FieldODEStateAndDerivative state0, T t) { + } + + public T g(FieldODEStateAndDerivative state) { + return state.getTime().subtract(12.0); + } + + public Action eventOccurred(FieldODEStateAndDerivative state, boolean increasing) { + return Action.STOP; + } + + public FieldODEState resetState(FieldODEStateAndDerivative state) { + return state; + } + + } + +} diff --git a/src/test/java/org/apache/commons/math3/ode/TestFieldProblem5.java b/src/test/java/org/apache/commons/math3/ode/TestFieldProblem5.java new file mode 100644 index 000000000..f84781ace --- /dev/null +++ b/src/test/java/org/apache/commons/math3/ode/TestFieldProblem5.java @@ -0,0 +1,41 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +package org.apache.commons.math3.ode; + +import org.apache.commons.math3.Field; +import org.apache.commons.math3.RealFieldElement; + +/** + * This class is used in the junit tests for the ODE integrators. + *

This is the same as problem 1 except integration is done + * backward in time

+ * @param the type of the field elements + */ +public class TestFieldProblem5> + extends TestFieldProblem1 { + + /** + * Simple constructor. + * @param field field to which elements belong + */ + public TestFieldProblem5(Field field) { + super(field); + setFinalConditions(getInitialTime().multiply(2).subtract(getFinalTime())); + } + +} diff --git a/src/test/java/org/apache/commons/math3/ode/TestFieldProblem6.java b/src/test/java/org/apache/commons/math3/ode/TestFieldProblem6.java new file mode 100644 index 000000000..cb4cd3eb6 --- /dev/null +++ b/src/test/java/org/apache/commons/math3/ode/TestFieldProblem6.java @@ -0,0 +1,84 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +package org.apache.commons.math3.ode; + +import org.apache.commons.math3.Field; +import org.apache.commons.math3.RealFieldElement; +import org.apache.commons.math3.util.MathArrays; + +/** + * This class is used in the junit tests for the ODE integrators. + + *

This specific problem is the following differential equation : + * 

+ *    y' = 3x^5 - y
+ *
+ * when the initial condition is y(0) = -360, the solution of this + * equation degenerates to a simple quintic polynomial function : + * 
+ *   y (t) = 3x^5 - 15x^4 + 60x^3 - 180x^2 + 360x - 360
+ *
+ *

+ + * @param the type of the field elements + */ +public class TestFieldProblem6> + extends TestFieldProblemAbstract { + + /** + * Simple constructor. + * @param field field to which elements belong + */ + public TestFieldProblem6(Field field) { + super(field); + setInitialConditions(convert(0.0), convert( new double[] { -360.0 })); + setFinalConditions(convert(1.0)); + setErrorScale(convert( new double[] { 1.0 })); + } + + @Override + public T[] doComputeDerivatives(T t, T[] y) { + + final T[] yDot = MathArrays.buildArray(getField(), getDimension()); + + // compute the derivatives + T t2 = t.multiply(t); + T t4 = t2.multiply(t2); + T t5 = t4.multiply(t); + for (int i = 0; i < getDimension(); ++i) { + yDot[i] = t5.multiply(3).subtract(y[i]); + } + + return yDot; + + } + + @Override + public T[] computeTheoreticalState(T t) { + + final T[] y = MathArrays.buildArray(getField(), getDimension()); + + for (int i = 0; i < getDimension(); ++i) { + y[i] = t.multiply(3).subtract(15).multiply(t).add(60).multiply(t).subtract(180).multiply(t).add(360).multiply(t).subtract(360); + } + + return y; + + } + +} diff --git a/src/test/java/org/apache/commons/math3/ode/TestFieldProblemAbstract.java b/src/test/java/org/apache/commons/math3/ode/TestFieldProblemAbstract.java new file mode 100644 index 000000000..09a10d810 --- /dev/null +++ b/src/test/java/org/apache/commons/math3/ode/TestFieldProblemAbstract.java @@ -0,0 +1,209 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +package org.apache.commons.math3.ode; + +import java.lang.reflect.Array; + +import org.apache.commons.math3.Field; +import org.apache.commons.math3.RealFieldElement; +import org.apache.commons.math3.ode.events.FieldEventHandler; +import org.apache.commons.math3.util.MathArrays; + +/** + * This class is used as the base class of the problems that are + * integrated during the junit tests for the ODE integrators. + * @param the type of the field elements + */ +public abstract class TestFieldProblemAbstract> + implements FieldFirstOrderDifferentialEquations { + + /** Field to which elements belong. */ + private Field field; + + /** Dimension of the problem. */ + private int n; + + /** Number of functions calls. */ + private int calls; + + /** Initial time */ + private T t0; + + /** Initial state */ + private T[] y0; + + /** Final time */ + private T t1; + + /** Error scale */ + private T[] errorScale; + + /** + * Simple constructor. + * @param field field to which elements belong + */ + protected TestFieldProblemAbstract(Field field) { + this.field = field; + n = 0; + calls = 0; + t0 = field.getZero(); + y0 = null; + t1 = field.getZero(); + errorScale = null; + } + + /** + * Set the initial conditions + * @param t0 initial time + * @param y0 initial state vector + */ + protected void setInitialConditions(T t0, T[] y0) { + calls = 0; + n = y0.length; + this.t0 = t0; + this.y0 = y0.clone(); + } + + /** + * Set the final conditions. + * @param t1 final time + */ + protected void setFinalConditions(T t1) { + this.t1 = t1; + } + + /** + * Set the error scale + * @param errorScale error scale + */ + protected void setErrorScale(T[] errorScale) { + this.errorScale = errorScale.clone(); + } + + /** get the filed to which elements belong. + * @return field to which elements belong + */ + public Field getField() { + return field; + } + + /** Get the problem dimension. + * @return problem dimension + */ + public int getDimension() { + return n; + } + + /** + * Get the initial time. + * @return initial time + */ + public T getInitialTime() { + return t0; + } + + /** + * Get the initial state vector. + * @return initial state vector + */ + public T[] getInitialState() { + return y0; + } + + /** + * Get the final time. + * @return final time + */ + public T getFinalTime() { + return t1; + } + + /** + * Get the error scale. + * @return error scale + */ + public T[] getErrorScale() { + return errorScale; + } + + /** + * Get the events handlers. + * @return events handlers */ + public FieldEventHandler[] getEventsHandlers() { + @SuppressWarnings("unchecked") + final FieldEventHandler[] empty = + (FieldEventHandler[]) Array.newInstance(FieldEventHandler.class, 0); + return empty; + } + + /** + * Get the theoretical events times. + * @return theoretical events times + */ + public T[] getTheoreticalEventsTimes() { + return MathArrays.buildArray(field, 0); + } + + /** + * Get the number of calls. + * @return nuber of calls + */ + public int getCalls() { + return calls; + } + + /** {@inheritDoc} */ + public void init(T t0, T[] y0, T t) { + } + + /** {@inheritDoc} */ + public T[] computeDerivatives(T t, T[] y) { + ++calls; + return doComputeDerivatives(t, y); + } + + abstract public T[] doComputeDerivatives(T t, T[] y); + + /** + * Compute the theoretical state at the specified time. + * @param t time at which the state is required + * @return state vector at time t + */ + abstract public T[] computeTheoreticalState(T t); + + /** Convert a double. + * @param d double to convert + * @return converted double + */ + protected T convert(double d) { + return field.getZero().add(d); + } + + /** Convert a one dimension array. + * @param elements array elements + * @return converted array + */ + protected T[] convert(double ... elements) { + T[] array = MathArrays.buildArray(field, elements.length); + for (int i = 0; i < elements.length; ++i) { + array[i] = convert(elements[i]); + } + return array; + } + +} diff --git a/src/test/java/org/apache/commons/math3/ode/TestFieldProblemHandler.java b/src/test/java/org/apache/commons/math3/ode/TestFieldProblemHandler.java new file mode 100644 index 000000000..d4e20e606 --- /dev/null +++ b/src/test/java/org/apache/commons/math3/ode/TestFieldProblemHandler.java @@ -0,0 +1,154 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +package org.apache.commons.math3.ode; + +import org.apache.commons.math3.RealFieldElement; +import org.apache.commons.math3.exception.MaxCountExceededException; +import org.apache.commons.math3.ode.sampling.FieldStepHandler; +import org.apache.commons.math3.ode.sampling.FieldStepInterpolator; +import org.apache.commons.math3.util.MathUtils; + +/** + * This class is used to handle steps for the test problems + * integrated during the junit tests for the ODE integrators. + * @param the type of the field elements + */ +public class TestFieldProblemHandler> + implements FieldStepHandler { + + /** Associated problem. */ + private TestFieldProblemAbstract problem; + + /** Maximal errors encountered during the integration. */ + private T maxValueError; + private T maxTimeError; + + /** Error at the end of the integration. */ + private T lastError; + + /** Time at the end of integration. */ + private T lastTime; + + /** ODE solver used. */ + private FieldFirstOrderIntegrator integrator; + + /** Expected start for step. */ + private T expectedStepStart; + + /** + * Simple constructor. + * @param problem problem for which steps should be handled + * @param integrator ODE solver used + */ + public TestFieldProblemHandler(TestFieldProblemAbstract problem, FieldFirstOrderIntegrator integrator) { + this.problem = problem; + this.integrator = integrator; + maxValueError = problem.getField().getZero(); + maxTimeError = problem.getField().getZero(); + lastError = problem.getField().getZero(); + expectedStepStart = null; + } + + public void init(FieldODEStateAndDerivative state0, T t) { + maxValueError = problem.getField().getZero(); + maxTimeError = problem.getField().getZero(); + lastError = problem.getField().getZero(); + expectedStepStart = null; + } + + public void handleStep(FieldStepInterpolator interpolator, boolean isLast) throws MaxCountExceededException { + + T start = integrator.getCurrentStepStart().getTime(); + if (start.subtract(problem.getInitialTime()).divide(integrator.getCurrentSignedStepsize()).abs().getReal() > 0.001) { + // multistep integrators do not handle the first steps themselves + // so we have to make sure the integrator we look at has really started its work + if (expectedStepStart != null) { + // the step should either start at the end of the integrator step + // or at an event if the step is split into several substeps + T stepError = MathUtils.max(maxTimeError, start.subtract(expectedStepStart).abs()); + for (T eventTime : problem.getTheoreticalEventsTimes()) { + stepError = MathUtils.min(stepError, start.subtract(eventTime).abs()); + } + maxTimeError = MathUtils.max(maxTimeError, stepError); + } + expectedStepStart = start.add(integrator.getCurrentSignedStepsize()); + } + + T pT = interpolator.getPreviousState().getTime(); + T cT = interpolator.getCurrentState().getTime(); + T[] errorScale = problem.getErrorScale(); + + // store the error at the last step + if (isLast) { + T[] interpolatedY = interpolator.getInterpolatedState(cT).getState(); + T[] theoreticalY = problem.computeTheoreticalState(cT); + for (int i = 0; i < interpolatedY.length; ++i) { + T error = interpolatedY[i].subtract(theoreticalY[i]).abs(); + lastError = MathUtils.max(error, lastError); + } + lastTime = cT; + } + + // walk through the step + for (int k = 0; k <= 20; ++k) { + + T time = pT.add(cT.subtract(pT).multiply(k).divide(20)); + T[] interpolatedY = interpolator.getInterpolatedState(time).getState(); + T[] theoreticalY = problem.computeTheoreticalState(time); + + // update the errors + for (int i = 0; i < interpolatedY.length; ++i) { + T error = errorScale[i].multiply(interpolatedY[i].subtract(theoreticalY[i]).abs()); + maxValueError = MathUtils.max(error, maxValueError); + } + } + } + + /** + * Get the maximal value error encountered during integration. + * @return maximal value error + */ + public T getMaximalValueError() { + return maxValueError; + } + + /** + * Get the maximal time error encountered during integration. + * @return maximal time error + */ + public T getMaximalTimeError() { + return maxTimeError; + } + + /** + * Get the error at the end of the integration. + * @return error at the end of the integration + */ + public T getLastError() { + return lastError; + } + + /** + * Get the time at the end of the integration. + * @return time at the end of the integration. + */ + public T getLastTime() { + return lastTime; + } + +}