はじめに

Javaで数値解析をしてみたかったのでapache commons mathライブラリの機能を使って実装してみた。
実装したのは高校物理でよく扱う斜方投射とした。

プログラム

Test.java

import java.util.List;

import org.apache.commons.math4.legacy.ode.FirstOrderDifferentialEquations;
import org.apache.commons.math4.legacy.ode.FirstOrderIntegrator;
import org.apache.commons.math4.legacy.ode.nonstiff.ClassicalRungeKuttaIntegrator;
import org.apache.commons.math4.legacy.exception.DimensionMismatchException;
import org.apache.commons.math4.legacy.exception.MaxCountExceededException;

public class Test {
    public static void main(String[] args) {
        // 運動方程式(1階微分方程式)
        FirstOrderDifferentialEquations ode = new FirstOrderDifferentialEquations() {
            //yの次元
            @Override
            public int getDimension() {
                return 4;
            }

            //1階微分方程式
            @Override
            public void computeDerivatives(double t, double[] y, double[] yDot)
                    throws MaxCountExceededException, DimensionMismatchException {
                double vx = y[0];
                double vy = y[1];

                double g = 9.8065; // 重力加速度
                yDot[0] = 0; // ax
                yDot[1] = -g; // ay
                yDot[2] = vx; // vx
                yDot[3] = vy; // vy
            }
        };
        
        //記録用のstpeHandlerを作成
        MyStepHandler stepHandler = new MyStepHandler();

        // 積分器Integratorインスタンス作成
        double step = 0.1 //0.1秒幅で計算
        FirstOrderIntegrator integrator = new ClassicalRungeKuttaIntegrator(step);
        integrator.addStepHandler(stepHandler);

        // 初期条件
        double startTime = 0;
        double endTime = 2;
        double v0 = 10.0;
        double theta = Math.toRadians(45.);
        double v0x = v0 * Math.cos(theta);
        double v0y = v0 * Math.sin(theta);
        double p0x = 0;
        double p0y = 0;
        double[] initialValues = new double[] { v0x, v0y, p0x, p0y };
        double[] endValues = new double[4];

        // 計算
        integrator.integrate(ode, startTime, initialValues, endTime, endValues);

        // 結果取得
        List<double[]> valueList = stepHandler.getValueList();
        List<Double> timeList = stepHandler.getTimeList();

        // 表示
        for (int i = 0; i < timeList.size(); i++) {
            double t = timeList.get(i);
            double[] y = valueList.get(i);
            double vx = y[0];
            double vy = y[1];
            double px = y[2];
            double py = y[3];
            System.out.println(String.format("%f\t%f\t%f\t%f\t%f\t", t, vx, vy, px, py));
        }
    }
}

MyStepHandler.java

import java.util.ArrayList;
import java.util.List;

import org.apache.commons.math4.legacy.ode.sampling.StepHandler;
import org.apache.commons.math4.legacy.ode.sampling.StepInterpolator;
import org.apache.commons.math4.legacy.exception.MaxCountExceededException;

public class MyStepHandler implements StepHandler {
    private List<double[]> valueList = new ArrayList<double[]>();
    private List<Double> timeList = new ArrayList<Double>();

    //計算開始時に呼び出されるメソッド
    @Override
    public void init(double t0, double[] y0, double t) {
        valueList.clear();
        timeList.clear();

        valueList.add(y0.clone());
        timeList.add(t0);
    }

    //各時刻の計算が終わるたびに呼び出されるメソッド
    @Override
    public void handleStep(StepInterpolator interpolator, boolean isLast) throws MaxCountExceededException {
        double[] y = interpolator.getInterpolatedState();
        double t = interpolator.getCurrentTime();

        valueList.add(y.clone());// クローンしないと値が上書きされる。
        timeList.add(t);
    }

    public List<double[]> getValueList() {
        return valueList;
    }

    public List<Double> getTimeList() {
        return timeList;
    }
};

結果

0.000000        7.071068        7.071068        0.000000        0.000000
0.100000        7.071068        6.090418        0.707107        0.658074
0.200000        7.071068        5.109768        1.414214        1.218084
0.300000        7.071068        4.129118        2.121320        1.680028
0.400000        7.071068        3.148468        2.828427        2.043907
0.500000        7.071068        2.167818        3.535534        2.309721
0.600000        7.071068        1.187168        4.242641        2.477471
0.700000        7.071068        0.206518        4.949747        2.547155
0.800000        7.071068        -0.774132       5.656854        2.518774
0.900000        7.071068        -1.754782       6.363961        2.392329
1.000000        7.071068        -2.735432       7.071068        2.167818
1.100000        7.071068        -3.716082       7.778175        1.845242
1.200000        7.071068        -4.696732       8.485281        1.424601
1.300000        7.071068        -5.677382       9.192388        0.905896
1.400000        7.071068        -6.658032       9.899495        0.289125
1.500000        7.071068        -7.638682       10.606602       -0.425711
1.600000        7.071068        -8.619332       11.313708       -1.238612
1.700000        7.071068        -9.599982       12.020815       -2.149577
1.800000        7.071068        -10.580632      12.727922       -3.158608
1.900000        7.071068        -11.561282      13.435029       -4.265704
2.000000        7.071068        -12.541932      14.142136       -5.470864

↓グラフにするとこんな感じ。

所感

ルンゲ・クッタ法がライブラリとしてあるのはありがたい。
StepHandlerを作成して、FirstOrderIntegratorにaddStepHandlerしないと各時刻ごとの計算結果を取得できないのがめんどくさい。

参考

公式ドキュメント
https://commons.apache.org/proper/commons-math/userguide/ode.html
ルンゲ・クッタ法について
https://ja.wikipedia.org/wiki/%E3%83%AB%E3%83%B3%E3%82%B2%EF%BC%9D%E3%82%AF%E3%83%83%E3%82%BF%E6%B3%95