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(ソイッシュ)を飲んだ。 »
「学問・資格」カテゴリの記事
- Wordleを100回やってみて正解率の分布を見た。1回目で当たる確率をpとして、2回目からはk倍ずつ確率が上がっていく(kp,k²p,k³p…)とすると、私の結果はp=0.0089, k=2.4になった。(2022.05.27)
- 高周波回路シミュレータQucsStudioを使ってみる(その5)マルチポート(4ポートだとs4pみたいな)のSパラメータデータをTouchstoneファイルに出力しようとするとExport to SnPでは2ポートになってしまう。Project→Import Dataからポート数を手で入れないと。(2022.05.24)
- 高周波(RF・マイクロ波・ミリ波・5G)関連ニュース2022年5月19日 IEEE Microwave Magazineで磁性体を使わない非可逆デバイス、Microwave Journalは5Gスモールセル、THz、Movanoがミリ波血圧&血糖値計発表、三菱の3Dプリンタ衛星アンテナ、など。(2022.05.19)
- シン・ウルトラマンをIMAXで観てきた。面白かった!物理学者の有岡君がラグランジアン,AdS,プランクブレーンとかしゃべってる!Randall–Sundrumも。マグカップは標準理論+重力で、物理監修は予想通り橋本幸士さんでした!ブラックホールの構造もタイムリー!(2022.05.14)
- 個人の正解確率pとして多数決をとったときの正解確率が、pが0.5に近づくとどのくらいの人数で0.9を超えるか計算してみた。p=0.50005だと日本の人口くらいの多数決。(神とさざなみの密室読んで正当性確率が気になったので)(2022.05.10)
« 拡散方程式をDOP853(ドルマンプリンス-ルンゲクッタ8次)で計算(Excel VBA) | トップページ | SOYSH(ソイッシュ)を飲んだ。 »
コメント