Как решить ОДУ с помощью Java?

Я пытаюсь решить ODE с помощью Java, и до сих пор я пробовал две разные библиотеки. Больше всего я доверяю Apache Commons Math, однако даже для простых задач я не вижу правильного решения.

Когда я решаю свою систему в Mathematica, я получаю следующее:

правильно

тогда как если я решу это с помощью решателя Дорманда-Принца 8 (5,3) в Apache Commons Math, я получу это:

неверно

Согласно теории, стоящей за этим, круг должен быть замкнут, как на первом рисунке. Как я могу получить это решение на 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]);
    }

}

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

}

person C. E.    schedule 12.03.2014    source источник


Ответы (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);

и посмотрите, приблизитесь ли вы к идеальному кругу.

Кстати: вычисление 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);
person Marco13    schedule 12.03.2014
comment
Вот почему я использую только python, потому что в python никогда не возникает путаницы с целочисленным делением... о, подождите. - person uhoh; 06.09.2017