mirror of
https://github.com/apache/commons-math.git
synced 2025-02-06 01:59:13 +00:00
Prepared tests for field-based ODE.
This commit is contained in:
Luc Maisonobe
parent
4b698fbf79
commit
7954482001
@ -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.
|
||||
|
||||
* <p>This specific problem is the following differential equation :
|
||||
* <pre>
|
||||
* y' = -y
|
||||
* </pre>
|
||||
* the solution of this equation is a simple exponential function :
|
||||
* <pre>
|
||||
* y (t) = y (t0) exp (t0-t)
|
||||
* </pre>
|
||||
* </p>
|
||||
|
||||
* @param <T> the type of the field elements
|
||||
*/
|
||||
public class TestFieldProblem1<T extends RealFieldElement<T>>
|
||||
extends TestFieldProblemAbstract<T> {
|
||||
|
||||
/**
|
||||
* Simple constructor.
|
||||
* @param field field to which elements belong
|
||||
*/
|
||||
public TestFieldProblem1(Field<T> 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;
|
||||
}
|
||||
|
||||
}
|
||||
@ -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.
|
||||
|
||||
* <p>This specific problem is the following differential equation :
|
||||
* <pre>
|
||||
* y' = t^3 - t y
|
||||
* </pre>
|
||||
* with the initial condition y (0) = 0. The solution of this equation
|
||||
* is the following function :
|
||||
* <pre>
|
||||
* y (t) = t^2 + 2 (exp (- t^2 / 2) - 1)
|
||||
* </pre>
|
||||
* </p>
|
||||
|
||||
* @param <T> the type of the field elements
|
||||
*/
|
||||
public class TestFieldProblem2<T extends RealFieldElement<T>>
|
||||
extends TestFieldProblemAbstract<T> {
|
||||
|
||||
/**
|
||||
* Simple constructor.
|
||||
* @param field field to which elements belong
|
||||
*/
|
||||
public TestFieldProblem2(Field<T> 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;
|
||||
}
|
||||
|
||||
}
|
||||
@ -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.
|
||||
|
||||
* <p>This specific problem is the following differential equation :
|
||||
* <pre>
|
||||
* 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
|
||||
* </pre>
|
||||
* This is a two-body problem in the plane which can be solved by
|
||||
* Kepler's equation
|
||||
* <pre>
|
||||
* y1 (t) = ...
|
||||
* </pre>
|
||||
* </p>
|
||||
|
||||
* @param <T> the type of the field elements
|
||||
*/
|
||||
public class TestFieldProblem3<T extends RealFieldElement<T>>
|
||||
extends TestFieldProblemAbstract<T> {
|
||||
|
||||
/** Eccentricity */
|
||||
T e;
|
||||
|
||||
/**
|
||||
* Simple constructor.
|
||||
* @param field field to which elements belong
|
||||
* @param e eccentricity
|
||||
*/
|
||||
public TestFieldProblem3(Field<T> 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<T> 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;
|
||||
|
||||
}
|
||||
|
||||
}
|
||||
@ -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.
|
||||
|
||||
* <p>This specific problem is the following differential equation :
|
||||
* <pre>
|
||||
* x'' = -x
|
||||
* </pre>
|
||||
* And when x decreases down to 0, the state should be changed as follows :
|
||||
* <pre>
|
||||
* x' -> -x'
|
||||
* </pre>
|
||||
* The theoretical solution of this problem is x = |sin(t+a)|
|
||||
* </p>
|
||||
|
||||
* @param <T> the type of the field elements
|
||||
*/
|
||||
public class TestFieldProblem4<T extends RealFieldElement<T>>
|
||||
extends TestFieldProblemAbstract<T> {
|
||||
|
||||
/** Time offset. */
|
||||
private T a;
|
||||
|
||||
/** Simple constructor.
|
||||
* @param field field to which elements belong
|
||||
*/
|
||||
public TestFieldProblem4(Field<T> 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<T>[] getEventsHandlers() {
|
||||
@SuppressWarnings("unchecked")
|
||||
FieldEventHandler<T>[] handlers =
|
||||
(FieldEventHandler<T>[]) Array.newInstance(FieldEventHandler.class, 2);
|
||||
handlers[0] = new Bounce<T>();
|
||||
handlers[1] = new Stop<T>();
|
||||
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<T extends RealFieldElement<T>> implements FieldEventHandler<T> {
|
||||
|
||||
private int sign;
|
||||
|
||||
public Bounce() {
|
||||
sign = +1;
|
||||
}
|
||||
|
||||
public void init(FieldODEStateAndDerivative<T> state0, T t) {
|
||||
}
|
||||
|
||||
public T g(FieldODEStateAndDerivative<T> state) {
|
||||
return state.getState()[0].multiply(sign);
|
||||
}
|
||||
|
||||
public Action eventOccurred(FieldODEStateAndDerivative<T> state, boolean increasing) {
|
||||
// this sign change is needed because the state will be reset soon
|
||||
sign = -sign;
|
||||
return Action.RESET_STATE;
|
||||
}
|
||||
|
||||
public FieldODEState<T> resetState(FieldODEStateAndDerivative<T> state) {
|
||||
T[] y = state.getState();
|
||||
y[0] = y[0].negate();
|
||||
y[1] = y[1].negate();
|
||||
return new FieldODEState<T>(state.getTime(), y);
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
private static class Stop<T extends RealFieldElement<T>> implements FieldEventHandler<T> {
|
||||
|
||||
public Stop() {
|
||||
}
|
||||
|
||||
public void init(FieldODEStateAndDerivative<T> state0, T t) {
|
||||
}
|
||||
|
||||
public T g(FieldODEStateAndDerivative<T> state) {
|
||||
return state.getTime().subtract(12.0);
|
||||
}
|
||||
|
||||
public Action eventOccurred(FieldODEStateAndDerivative<T> state, boolean increasing) {
|
||||
return Action.STOP;
|
||||
}
|
||||
|
||||
public FieldODEState<T> resetState(FieldODEStateAndDerivative<T> state) {
|
||||
return state;
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
}
|
||||
@ -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.
|
||||
* <p>This is the same as problem 1 except integration is done
|
||||
* backward in time</p>
|
||||
* @param <T> the type of the field elements
|
||||
*/
|
||||
public class TestFieldProblem5<T extends RealFieldElement<T>>
|
||||
extends TestFieldProblem1<T> {
|
||||
|
||||
/**
|
||||
* Simple constructor.
|
||||
* @param field field to which elements belong
|
||||
*/
|
||||
public TestFieldProblem5(Field<T> field) {
|
||||
super(field);
|
||||
setFinalConditions(getInitialTime().multiply(2).subtract(getFinalTime()));
|
||||
}
|
||||
|
||||
}
|
||||
@ -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.
|
||||
|
||||
* <p>This specific problem is the following differential equation :
|
||||
* <pre>
|
||||
* y' = 3x^5 - y
|
||||
* </pre>
|
||||
* when the initial condition is y(0) = -360, the solution of this
|
||||
* equation degenerates to a simple quintic polynomial function :
|
||||
* <pre>
|
||||
* y (t) = 3x^5 - 15x^4 + 60x^3 - 180x^2 + 360x - 360
|
||||
* </pre>
|
||||
* </p>
|
||||
|
||||
* @param <T> the type of the field elements
|
||||
*/
|
||||
public class TestFieldProblem6<T extends RealFieldElement<T>>
|
||||
extends TestFieldProblemAbstract<T> {
|
||||
|
||||
/**
|
||||
* Simple constructor.
|
||||
* @param field field to which elements belong
|
||||
*/
|
||||
public TestFieldProblem6(Field<T> 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;
|
||||
|
||||
}
|
||||
|
||||
}
|
||||
@ -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 <T> the type of the field elements
|
||||
*/
|
||||
public abstract class TestFieldProblemAbstract<T extends RealFieldElement<T>>
|
||||
implements FieldFirstOrderDifferentialEquations<T> {
|
||||
|
||||
/** Field to which elements belong. */
|
||||
private Field<T> 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<T> 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<T> 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<T>[] getEventsHandlers() {
|
||||
@SuppressWarnings("unchecked")
|
||||
final FieldEventHandler<T>[] empty =
|
||||
(FieldEventHandler<T>[]) 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;
|
||||
}
|
||||
|
||||
}
|
||||
@ -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 <T> the type of the field elements
|
||||
*/
|
||||
public class TestFieldProblemHandler<T extends RealFieldElement<T>>
|
||||
implements FieldStepHandler<T> {
|
||||
|
||||
/** Associated problem. */
|
||||
private TestFieldProblemAbstract<T> 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<T> 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<T> problem, FieldFirstOrderIntegrator<T> 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<T> state0, T t) {
|
||||
maxValueError = problem.getField().getZero();
|
||||
maxTimeError = problem.getField().getZero();
|
||||
lastError = problem.getField().getZero();
|
||||
expectedStepStart = null;
|
||||
}
|
||||
|
||||
public void handleStep(FieldStepInterpolator<T> 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;
|
||||
}
|
||||
|
||||
}
|
||||
Loading…
x
Reference in New Issue
Block a user