Javaの数値計算ライブラリApache Commons Mathを使う(5) 常微分方程式の数値解法、アダプティブな刻み幅の8次のルンゲクッタ法、Dormand-Princeを使ってローレンツ方程式を計算、 JFreeChartでプロット。

今回はこちらの例題から。

 Visual C# (C_sharp)の数学ライブラリ Math.NET Numericsを使う(5) 常微分方程式の数値解法、4段4次のルンゲクッタ法がRungeKutta.FourthOrderの一文でできる。ローレンツ方程式を例としてやってみる。

ローレンツ方程式はこんな形。

\begin{align}
&\frac{dx}{dt} = \sigma (y - x) \\\
&\frac{dy}{dt} = -xz + rx -y \\\
&\frac{dz}{dt} = xy - bz \\\
\end{align}

Appache Commons Mathはいろいろな常微分方程式用のソルバーがある。

Fixed Step Integrators
Name	Order
Euler	1
Midpoint	2
Classical Runge-Kutta	4
Gill	4
3/8	4
Luther	6

Adaptive Stepsize Integrators
Name	Integration Order	Error Estimation Order
Higham and Hall	5	4
Dormand-Prince 5(4)	5	4
Dormand-Prince 8(5,3)	8	5 and 3
Gragg-Bulirsch-Stoer	variable (up to 18 by default)	variable
Adams-Bashforth	variable	variable
Adams-Moulton	variable	variable

ここはうちのブログでよく使っている8次のDormand Princeにしておこう。

コードはこんな感じで。JFreeChatには罠があって、そのままプロットするとX軸の小さい順にプロットしてしまう。XYSeriesの引数にfalseを付けるとそのままのデータでプロットしてくれる。

import java.awt.BorderLayout;

import java.util.ArrayList;

import javax.swing.JFrame;

import org.apache.commons.math3.ode.FirstOrderDifferentialEquations;

import org.apache.commons.math3.ode.FirstOrderIntegrator;

import org.apache.commons.math3.ode.nonstiff.DormandPrince853Integrator;

import org.apache.commons.math3.ode.sampling.FixedStepHandler;

import org.apache.commons.math3.ode.sampling.StepHandler;

import org.apache.commons.math3.ode.sampling.StepInterpolator;

import org.apache.commons.math3.ode.sampling.StepNormalizer;

import org.jfree.chart.ChartFactory;

import org.jfree.chart.ChartPanel;

import org.jfree.chart.JFreeChart;

import org.jfree.chart.plot.PlotOrientation;

import org.jfree.chart.plot.XYPlot;

import org.jfree.chart.renderer.xy.XYLineAndShapeRenderer;

import org.jfree.data.xy.XYSeries;

import org.jfree.data.xy.XYSeriesCollection;

public class LorenzEquation extends JFrame {

    private static final long serialVersionUID = 1L;

    public static void main(String[] args) {

        LorenzEquation frame = new LorenzEquation();

        frame.setDefaultCloseOperation(JFrame.EXIT_ON_CLOSE);

        frame.setBounds(10, 10, 640, 480);

        frame.setTitle("Lorenz Equation");

        frame.setVisible(true);

    }

    public LorenzEquation() {

         JFreeChart chart =

                  ChartFactory.createXYLineChart("Lorenz Equation by Dormand-Prince 8(5,3)",

                                                 "X",

                                                 "Y, Z",

                                                 createData(),

                                                 PlotOrientation.VERTICAL,

                                                 true,

                                                 false,

                                                 false);

                XYPlot plot = chart.getXYPlot();

                XYLineAndShapeRenderer renderer =new XYLineAndShapeRenderer();

                renderer.setSeriesShapesVisible(0, false);

                renderer.setSeriesShapesVisible(1, false);

                plot.setRenderer(renderer);

                ChartPanel cpanel = new ChartPanel(chart);

                getContentPane().add(cpanel, BorderLayout.CENTER);

    }

    private static class ODE implements FirstOrderDifferentialEquations {

        private double[] param;

        public ODE(double[] param) {

            this.param = param;

        }

        public int getDimension() {

            return 3;

        }

        public void computeDerivatives(double t, double[] y, double[] yDot) {

            yDot[0] = param[0] * (y[1] - y[0]);

            yDot[1] = param[1] * y[0] - y[1] - y[0] * y[2];

            yDot[2] = y[0] * y[1] - param[2] * y[2];

        }

    }

    private XYSeriesCollection createData(){

        ArrayList<Double[]> ylist = new ArrayList<>();

        ArrayList<Double> tlist = new ArrayList<>();

        double tstep = 0.01;

        double tmax = 200.0;

        FirstOrderIntegrator dp853 = new DormandPrince853Integrator(1.0e-8, tstep, 1.0e-10, 1.0e-10);

        FirstOrderDifferentialEquations ode = new ODE(new double[] { 10.0, 28.0, 8.0 / 3.0} );

        double[] y = new double[] { 1.0, 1.0, 1.0 }; // initial state

        StepHandler stepHandler = new StepHandler() {

            public void init(double t0, double[] y0, double t) {

            }

            public void handleStep(StepInterpolator interpolator, boolean isLast) {

                double   t = interpolator.getCurrentTime();

                double[] y = interpolator.getInterpolatedState();

                ylist.add(new Double[] {y[0], y[1], y[2]});

                tlist.add(t);

            }

        };

        dp853.addStepHandler(stepHandler);

        dp853.integrate(ode, 0.0, y, tmax, y);

        XYSeriesCollection data = new XYSeriesCollection();

        XYSeries series1 = new XYSeries("Y", false);

        XYSeries series2 = new XYSeries("Z", false);

        for (int i = 0 ; i < ylist.size(); i++){

          series1.add(ylist.get(i)[0], ylist.get(i)[1]);

          series2.add(ylist.get(i)[0], ylist.get(i)[2]);

        }

        data.addSeries(series1);

        data.addSeries(series2);

        return data;

      }

}

結果はこちら。