How to solve ODEs with Java?

Question

I'm trying to solve an ODE with Java and so far I have tried two different libraries. The one I trust the most is Apache Commons Math however even for simple problems I don't seem to get the correct solution.

When I solve my system in Mathematica I get this:

whereas if I solve it with the Dormand-Prince 8(5,3) solver in Apache Commons Math I get this:

According to the theory behind this the circle should be closed like in the first picture. How can I get this solution in Java?

Here is the Java code I currently use (I have tested different settings for step size etc.):

program.java:

import org.apache.commons.math3.ode.*;
import org.apache.commons.math3.ode.nonstiff.DormandPrince853Integrator;
import org.apache.commons.math3.ode.sampling.StepHandler;
import org.apache.commons.math3.ode.sampling.StepInterpolator;

import java.io.File;
import java.io.PrintWriter;
import java.util.ArrayList;

import static java.lang.Math.*;

public class program {

    public static void main(String[] args) {

        FirstOrderIntegrator dp853 = new DormandPrince853Integrator(1.0E-8, 10E-5, 1.0E-20, 1.0E-20);
        FirstOrderDifferentialEquations ode = new CometSun();

        double G = 1.48808E-34;
        double sunMass = 1.9891E30;

        double a = 1;
        double e = 0.1;
        double rp = a*(1-e);
        double v0 = sqrt(G*sunMass*(2/rp-1/a));
        double t = 2*PI*sqrt(pow(a,3)/(G*sunMass));

        double[] y = new double[6];

        y[0] = -rp;
        y[1] = 0;
        y[2] = 0;

        y[3] = 0;
        y[4] = v0;
        y[5] = 0;

        StepHandler stepHandler = new StepHandler() {

            ArrayList<String> steps = new ArrayList<String>();

            public void init(double t0, double[] y0, double t) {

            }

            public void handleStep(StepInterpolator interpolator, boolean isLast) {
                double   t = interpolator.getCurrentTime();
                double[] y = interpolator.getInterpolatedState();

                if( t > steps.size() )
                    steps.add(t + " " + y[0] + " " + y[1] + " " + y[2]);

                if(isLast) {
                    try{
                        PrintWriter writer = new PrintWriter(new File("results.txt"), "UTF-8");
                        for(String step: steps) {
                            writer.println(step);
                        }
                        writer.close();
                    } catch(Exception e) {};
                }
            }
        };
        dp853.addStepHandler(stepHandler);

        dp853.integrate(ode, 0.0, y, t, y);

        System.out.println(y[0]);
        System.out.println(y[1]);
    }

}

CometSun.java:

import org.apache.commons.math3.ode.FirstOrderDifferentialEquations;

import static java.lang.Math.*;

public class CometSun implements FirstOrderDifferentialEquations {

    double G = 1.48808E-34;
    double sunMass = 1.9891E30;

    public int getDimension() {
        return 6;
    }

    public void computeDerivatives(double t, double[] y, double[] yp) {
        double coeff = G*sunMass/pow(pow(y[0],2)+pow(y[1],2)+pow(y[2],2),3/2);

        yp[0] = y[3];
        yp[1] = y[4];
        yp[2] = y[5];

        yp[3] = -coeff*y[0];
        yp[4] = -coeff*y[1];
        yp[5] = -coeff*y[2];
    }

}

Answer 1

I can neither predict nor promise or verify that you'll receive exactly the same results as Mathematica: There are some very large and some very small numbers involved, and you might hit the limitations of the double precision.

But in this case, the main reason for the error is most likely a very simple and very subtle one:

double coeff = G*sunMass/pow(pow(y[0],2)+pow(y[1],2)+pow(y[2],2),3/2);

The last exponent is 1. You are performing an integer division there, and 3/2 yields 1. You could change this to

double coeff = G*sunMass/pow(pow(y[0],2)+pow(y[1],2)+pow(y[2],2),3.0/2.0);

and see whether you get closer to the perfect circle.

BTW: Computing x² as pow(x,2) is very inefficient. If efficiency may be an issue here (and I assume that it is), you should consider writing this as something like

double yy0 = y[0]*y[0];
double yy1 = y[1]*y[1];
double yy2 = y[2]*y[2];
double coeff = G*sunMass/pow(yy0+yy1+yy2,3.0/2.0);