ローレンツ方程式をExcel VBAでDormand-Prince8次で計算する（ソース付き）。

一周年ということで、今までのExcelで計算シリーズのソースでも公開していこうと思う。まずはRunge-Kutta8次であるDormand-Prince法でローレンツ方程式を計算したもの。

で、係数の表は以下の通り。これは埋め込み型の公式で、自動刻み幅調整もできるけどプログラムではパス（めんどくさい。遊びでやってるんだし）。

VBAのプログラムリストは以下の通り。適当にボタンを作って、それを押すと計算するようにしたらいいです。

「Lorenz-eq-Dormand-Prince8.txt」をダウンロード

パラメータはσ,r,bですが、変数としてはs,r,bbとしています。

関数f1,f2,f3を好きなものに変えるといろんな常微分方程式が計算できますよ。

計算結果をGnuplotで書くと、おなじみのグラフになります。

～　プログラムの中身　～

Option Explicit

Dim s As Double, r As Double, bb As Double

Private Sub CommandButton1_Click()
    Dim x As Double, y As Double, z As Double
    Dim xd As Double, yd As Double, zd As Double
   
    Dim kx(13) As Double
    Dim ky(13) As Double
    Dim kz(13) As Double
    Dim a(13) As Double
    Dim b(13, 13) As Double
    Dim c(13) As Double
       
    Dim t As Double, dt As Double
    Dim i As Long, j As Integer, k As Integer
    Dim pi As Double
    Application.ScreenUpdating = False
   
' パラメータ
    s = 10#
    r = 26#
    bb = 8# / 3#

'　初期値   
    x = 0.01
    y = 0.01
    z = 0
    t = 0
    dt = 1# / 100#
   
    a(1) = 0#
    a(2) = 2# / 27#
    a(3) = 1# / 9#
    a(4) = 1# / 6#
    a(5) = 5# / 12#
    a(6) = 0.5
    a(7) = 5# / 6#
    a(8) = 1# / 6#
    a(9) = 2# / 3#
    a(10) = 1# / 3#
    a(11) = 1#
    a(12) = 0#
    a(13) = 1
    For j = 1 To 13
        For k = 1 To 13
            b(j, k) = 0
        Next k
        kx(j) = 0#
        ky(j) = 0#
kz(j) = 0#
    Next j
   
    b(2, 1) = 2# / 27#
   
    b(3, 1) = 1# / 36
    b(3, 2) = 1# / 12#
   
    b(4, 1) = 1# / 24#
    b(4, 3) = 1# / 8#
   
    b(5, 1) = 5# / 12#
    b(5, 3) = -25# / 16#
    b(5, 4) = 25# / 16#
   
    b(6, 1) = 1# / 20#
    b(6, 4) = 1# / 4#
    b(6, 5) = 1# / 5#
   
    b(7, 1) = -25# / 108#
    b(7, 4) = 125# / 108#
    b(7, 5) = -65# / 27#
    b(7, 6) = 125# / 54#
   
    b(8, 1) = 31# / 300#
    b(8, 5) = 61# / 225#
    b(8, 6) = -2# / 9#
    b(8, 7) = 13# / 900#
   
    b(9, 1) = 2#
    b(9, 4) = -53# / 6#
    b(9, 5) = 704# / 45#
    b(9, 6) = -107# / 9#
    b(9, 7) = 67# / 90#
    b(9, 8) = 3#
   
    b(10, 1) = -91# / 108#
    b(10, 4) = 23# / 108#
    b(10, 5) = -976# / 135#
    b(10, 6) = 311# / 54#
    b(10, 7) = -19# / 60#
    b(10, 8) = 17# / 6#
    b(10, 9) = -1# / 12#
   
    b(11, 1) = 2383# / 4100#
    b(11, 4) = -341# / 164#
    b(11, 5) = 4496# / 1025#
    b(11, 6) = -301# / 82#
    b(11, 7) = 2133# / 4100#
    b(11, 8) = 45# / 82#
    b(11, 9) = 45# / 164#
    b(11, 10) = 18# / 41#
   
    b(12, 1) = 3# / 205#
    b(12, 6) = -6# / 41#
    b(12, 7) = -3# / 205#
    b(12, 8) = -3# / 41#
    b(12, 9) = 3# / 41#
    b(12, 10) = 6# / 41#
   
    b(13, 1) = -1777# / 4100#
    b(13, 4) = -341# / 164#
    b(13, 5) = 4496# / 1025#
    b(13, 6) = -289# / 82#
    b(13, 7) = 2193# / 4100#
    b(13, 8) = 51# / 82#
    b(13, 9) = 33# / 164#
    b(13, 10) = 12# / 41#
    b(13, 12) = 1#
   
    c(6) = 34# / 105#
    c(7) = 9# / 35#
    c(8) = 9# / 35#
    c(9) = 9# / 280#
    c(10) = 9# / 280#
    c(12) = 41# / 840#
    c(13) = 41# / 840#
   
    
    For i = 0 To 10000

        Worksheets("Sheet1").Cells(i + 1, 3) = t
        Worksheets("Sheet1").Cells(i + 1, 4) = x
        Worksheets("Sheet1").Cells(i + 1, 5) = y
        Worksheets("Sheet1").Cells(i + 1, 6) = z
       
       
        kx(1) = f1(t + dt, x, y, z) * dt
        ky(1) = f2(t + dt, x, y, z) * dt
        kz(1) = f3(t + dt, x, y, z) * dt
       
        For j = 2 To 13
            xd = x
            yd = y
            zd = z
            For k = 1 To j - 1
                xd = xd + b(j, k) * kx(k)
                yd = yd + b(j, k) * ky(k)
                zd = zd + b(j, k) * kz(k)
            Next k
            kx(j) = f1(t + a(j) * dt, xd, yd, zd) * dt
            ky(j) = f2(t + a(j) * dt, xd, yd, zd) * dt
            kz(j) = f3(t + a(j) * dt, xd, yd, zd) * dt
        Next j
       
       
        For j = 1 To 13
            x = x + c(j) * kx(j)
            y = y + c(j) * ky(j)
            z = z + c(j) * kz(j)
        Next j

        t = t + dt
       
    Next i
    Application.ScreenUpdating = True
End Sub

Function f1(t As Double, x As Double, y As Double, z As Double) As Double
    f1 = -s * (x - y)
End Function

Function f2(t As Double, x As Double, y As Double, z As Double) As Double
    f2 = -y - x * z + r * x
End Function

Function f3(t As Double, x As Double, y As Double, z As Double) As Double
    f3 = x * y - bb * z
End Function