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

今回はこの例題。

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

これはChatGPTに全面的にサポートしてもらってやっとできた…コードはこんな感じで。

(ns mini.Lorenz

   (:import [org.apache.commons.math3.ode FirstOrderDifferentialEquations]

            [org.apache.commons.math3.ode.nonstiff DormandPrince853Integrator]

            [org.apache.commons.math3.ode.sampling StepHandler StepInterpolator]

            [org.jfree.data.xy XYSeries XYSeriesCollection]

            [org.jfree.chart ChartFactory ChartPanel]

            [org.jfree.chart.plot PlotOrientation]

            [org.jfree.chart.renderer.xy XYLineAndShapeRenderer]

            [javax.swing JFrame]

            [java.awt BorderLayout]))

 ;; Lorenz ODE: dx/dt = σ(y-x), dy/dt = ρx - y - xz, dz/dt = xy - βz

 (defn make-lorenz [sigma rho beta]

   (proxy [FirstOrderDifferentialEquations] []

     (getDimension [] 3)

     (computeDerivatives [t y yDot]

       (aset yDot 0 (* sigma (- (aget y 1) (aget y 0))))

       (aset yDot 1 (- (* rho (aget y 0))

                       (aget y 1)

                       (* (aget y 0) (aget y 2))))

       (aset yDot 2 (- (* (aget y 0) (aget y 1))

                       (* beta (aget y 2)))))))

(defn integrate-lorenz []

  (let [ylist (atom [])

        tstep 0.01

        tmax 200.0

        integrator (DormandPrince853Integrator. 1.0e-8 tstep 1.0e-10 1.0e-10)

        ode (make-lorenz 10.0 28.0 (/ 8.0 3.0))

        y0 (double-array [1.0 1.0 1.0])

        step-handler

        (proxy [StepHandler] []

          (init [t0 y0 t])

          (handleStep [^StepInterpolator interpolator isLast]

            (let [t (.getCurrentTime interpolator)

                  y (.getInterpolatedState interpolator)]

              (swap! ylist conj [t (aget y 0) (aget y 1) (aget y 2)]))))]

    (.addStepHandler integrator step-handler)

    (.integrate integrator ode 0.0 y0 tmax y0)

    @ylist))

(defn create-dataset [points]

  (let [seriesY (XYSeries. "Y" false)

        seriesZ (XYSeries. "Z" false)]

    #_{:clj-kondo/ignore [:unused-binding]}

    (doseq [[t x y z] points]

      (.add seriesY x y)

      (.add seriesZ x z))

    (doto (XYSeriesCollection.)

      (.addSeries seriesY)

      (.addSeries seriesZ))))

(defn create-chart [dataset]

  (let [chart (ChartFactory/createXYLineChart

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

               "X"

               "Y, Z"

               dataset

               PlotOrientation/VERTICAL

               true false false)

        plot (.getXYPlot chart)

        renderer (XYLineAndShapeRenderer.)]

    (.setSeriesShapesVisible renderer 0 false)

    (.setSeriesShapesVisible renderer 1 false)

    (.setRenderer plot renderer)

    chart))

(defn -main []

  (let [points (integrate-lorenz)

        dataset (create-dataset points)

        chart (create-chart dataset)

        frame (JFrame. "Lorenz Equation")]

    (.setDefaultCloseOperation frame JFrame/EXIT_ON_CLOSE)

    (.setBounds frame 10 10 640 480)

    (doto frame

      (.add (ChartPanel. chart) BorderLayout/CENTER)

      (.setVisible true))))

(-main)

結果はこちら。