java - 如何用 Java 求解 ODE？

Question

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

当我在 Mathematica 中解决我的系统时，我得到了这个:

而如果我在 Apache Commons Math 中使用 Dormand-Prince 8(5,3) 求解器求解它，我会得到:

根据这背后的理论，圆圈应该像第一张图片中那样闭合。如何在 Java 中获得此解决方案？

这是我目前使用的 Java 代码(我测试了步长等的不同设置):

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

}

cometd .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

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

但在这种情况下，错误的主要原因很可能是一个非常简单且非常微妙的原因:

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

看看你是否离完美的圆圈更近了。

顺便说一句:将 x² 计算为 pow(x,2) 非常效率低下。如果效率在这里可能是一个问题(我认为是)，您应该考虑将其写成类似

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