src/test/java/org/apache/commons/math3/ode/TestFieldProblem3.java

/*
 * 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;

    }

}