« 拡散方程式をDOP853(ドルマンプリンス-ルンゲクッタ8次)で計算(Excel VBA) | トップページ | SOYSH(ソイッシュ)を飲んだ。 »

2011年5月19日 (木)

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

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

Kdv1dop853

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

Kdv2dop853

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

« 拡散方程式をDOP853(ドルマンプリンス-ルンゲクッタ8次)で計算(Excel VBA) | トップページ | SOYSH(ソイッシュ)を飲んだ。 »

学問・資格」カテゴリの記事

コメント

コメントを書く

(ウェブ上には掲載しません)

トラックバック

この記事のトラックバックURL:
http://app.cocolog-nifty.com/t/trackback/512682/51645763

この記事へのトラックバック一覧です: KdV方程式をDOP853(ルンゲクッタ8次)でExcel VBAで計算:

« 拡散方程式をDOP853(ドルマンプリンス-ルンゲクッタ8次)で計算(Excel VBA) | トップページ | SOYSH(ソイッシュ)を飲んだ。 »

2018年6月
          1 2
3 4 5 6 7 8 9
10 11 12 13 14 15 16
17 18 19 20 21 22 23
24 25 26 27 28 29 30
フォト
無料ブログはココログ