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ニュース 2025年7月13日 Pythonの高周波ライブラリscikit-rfがv1.8.0に、SamsungがGalaxy Z Fold7など発表→QualcommがSnapdragon 8 Eliteが使われていると発表、NGMNが基地局アンテナの推奨事項をまとめる、STMicroとMetalenzがメタサーフェス光学のライセンス締結(2025.07.14)
- Interface2025年8月号Pythonで体験!はじめての暗号を買った。上杉暗号からRSA、AES、DHなど、特に楕円曲線暗号についてはコードも実際に動かすところまで詳しくかかれていた。耐量子暗号や聞いたことなかったY-00暗号や関数型暗号も記載。(2025.07.10)
- 高周波・RFニュース 2025年7月8日 NordicとSercommのセルラーIoTモジュール、iFixitがFairphone 6を分解、スコアは10/10、RCR wireless newsのウェビナー2件(6GとIndustry4.0)、SEMCOが高耐圧C0G MLCCを車載急速充電に提案(2025.07.09)
- 高周波・RFニュース 2025年7月2日 5G Americasが6Gに向けセンシングと通信ホワイトペーパー発行、KYOCERA AVXが3dBハイブリッドカプラ発表、TDKが車載薄膜インダクタ発表、Nordicが1次電池向けPMIC発表、ローデ・シュワルツの6GとAI/ML解説記事(2025.07.02)
- 高周波・RFニュース 2025年6月30日 QualcommがAIを用いた6Rxアンテナ解説、Next G Allianceと日本のXGMFが5G,6Gで協力、5G Americasが25Q1で5G加入者増加と発表、TechInsigtsがHuawei Pura 80 Pro+分解、Qorvoが5-7GHzのWi-Fi 7 FEM発表(2025.06.30)
« 拡散方程式をDOP853(ドルマンプリンス-ルンゲクッタ8次)で計算(Excel VBA) | トップページ | SOYSH(ソイッシュ)を飲んだ。 »
コメント