KdV方程式をDOP853（ルンゲクッタ8次）でExcel VBAで計算

ドルマン・プリンス（Dormand-Prince)のDOP853で偏微分方程式を計算しようシリーズ。

KdV方程式

∂u/∂t　+u*∂u/∂x+δ^3 ∂^3u/∂x^3 = 0

をやってみよう。DOP853に与える関数のルーチンはこんな感じで。

Sub FCN(N As Long, X As Double, Y() As Double, F() As Double, RPAR() As Double, IPAR() As Long)

        Dim dx As Double, dt As Double, diff1 As Double, diff3 As Double, u As Double
        Dim p1 As Long, m1 As Long, p2 As Long, m2 As Long
        Dim I As Long
        dx = RPAR(1)
        dt = RPAR(2)
        u = RPAR(3)
       
        For I = 1 To N
            If I < N Then
                p1 = I + 1
            Else
                p1 = 1
            End If
            If I > 1 Then
                m1 = I - 1
            Else
                m1 = N
            End If
            If I < N - 1 Then
                p2 = I + 2
            ElseIf I = N - 1 Then
                p2 = 1
            Else
                p2 = 2
            End If
            If I > 2 Then
                m2 = I - 2
            ElseIf I = 2 Then
                m2 = N
            Else
                m2 = N - 1
            End If
            
            diff3 = (-Y(m2) + 2 * Y(m1) - 2 * Y(p1) + Y(p2)) / (2 * (dx ^ 3))
            diff1 = (Y(m2) - 8 * Y(m1) + 8 * Y(p1) - Y(p2)) / (12 * dx)
            F(I) = -Y(I) * diff1 - u * diff3
        Next I

End Sub

で計算結果はこちら。オリジナルのザブスキーとクルスカルのとパラメータはあわせてある。

さらに時間を進めたときの結果はこちら。

本当に何も考えなくても計算できる。便利。