如何用java解决ODEs？

Question

我正在尝试用Java解决一个ODE,到目前为止,我已经尝试了两个不同的库。我最信任的是Apache Commons Math，然而，即使是简单的问题，我似乎也得不到正确的解。

当我在Mathematica中求解系统时,我得到了这样的结果：

如果我用Apache Commons Math中的Dormand-Prince8(5,3）求解器求解,我会得到以下结果：

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]);
    }

}

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

我既不能预测，也不能保证或证实您会得到与Mathematica完全相同的结果：其中涉及一些非常大的数字和一些非常小的数字，您可能会遇到double精度的限制。

但在这种情况下，错误的主要原因极有可能是一个非常简单、非常微妙的原因：

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

最后一个指数是1。你在那里执行一个整数除法，3/2得到1。你可以把它改成

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

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);