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
で計算結果はこちら。オリジナルのザブスキーとクルスカルのとパラメータはあわせてある。
さらに時間を進めたときの結果はこちら。
本当に何も考えなくても計算できる。便利。
« 拡散方程式をDOP853(ドルマンプリンス-ルンゲクッタ8次)で計算(Excel VBA) | トップページ | SOYSH(ソイッシュ)を飲んだ。 »
「学問・資格」カテゴリの記事
- 高周波(RF・マイクロ波・ミリ波・5G)関連ニュース2021年2月16日 IEEE Microwave Magazineの特集はオールデジタルのRFID、Microwave JournalはEバンド ミリ波通信に衛星や気球を使う話、アメリカの半導体企業がバイデンに投資を迫る、(2021.02.17)
- カオスを生じる電気回路、Chua’s circuitをLTspiceで回路シミュレーションしてみる。(2021.02.19)
- Labyrinth Chaos(迷宮カオス)を生むThomas-Rössler方程式のパラメータbを色々変えて、Python+Scipyでルンゲクッタ8次のDOP853(Dormand&Prince)を使って計算してGIFアニメ(2021.02.16)
- フィッツヒュー・南雲 (FitzHugh-Nagumo) 方程式をPython+Scipyでルンゲクッタ8次のDOP853(Dormand Prince)で計算。(2021.02.23)
- 「水晶振動子の等価回路計算」をカシオの高精度計算サイトkeisan.casio.jpの自作式としてUP! インピーダンスの大きさと位相がグラフ化できる。(2021.02.12)
« 拡散方程式をDOP853(ドルマンプリンス-ルンゲクッタ8次)で計算(Excel VBA) | トップページ | SOYSH(ソイッシュ)を飲んだ。 »
コメント