如何使用Java解决常微分方程？

Question

如何使用Java解决常微分方程？

javamathapache-commonsode

6

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

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

而如果我使用Apache Commons Math中的Dormand-Prince 8(5,3)求解器进行求解，我会得到以下结果：

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

这是我当前使用的Java代码（我测试了不同的步长等设置）： 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];
    }

}

- C. E.

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Marco13 · Accepted Answer

我既不能预测，也不能保证或验证你会得到与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);

看看是否能让你更接近完美的圆。

顺便说一下：计算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);