ローレンツ方程式をExcel VBAでDormand-Prince8次で計算する(ソース付き)。
一周年ということで、今までのExcelで計算シリーズのソースでも公開していこうと思う。まずはRunge-Kutta8次であるDormand-Prince法でローレンツ方程式を計算したもの。
で、係数の表は以下の通り。これは埋め込み型の公式で、自動刻み幅調整もできるけどプログラムではパス(めんどくさい。遊びでやってるんだし)。
VBAのプログラムリストは以下の通り。適当にボタンを作って、それを押すと計算するようにしたらいいです。
「Lorenz-eq-Dormand-Prince8.txt」をダウンロード
パラメータはσ,r,bですが、変数としてはs,r,bbとしています。
関数f1,f2,f3を好きなものに変えるといろんな常微分方程式が計算できますよ。
計算結果をGnuplotで書くと、おなじみのグラフになります。
~ プログラムの中身 ~
Option Explicit
Dim s As Double, r As Double, bb As Double
Private Sub CommandButton1_Click()
Dim x As Double, y As Double, z As Double
Dim xd As Double, yd As Double, zd As Double
Dim kx(13) As Double
Dim ky(13) As Double
Dim kz(13) As Double
Dim a(13) As Double
Dim b(13, 13) As Double
Dim c(13) As Double
Dim t As Double, dt As Double
Dim i As Long, j As Integer, k As Integer
Dim pi As Double
Application.ScreenUpdating = False
' パラメータ
s = 10#
r = 26#
bb = 8# / 3#
' 初期値
x = 0.01
y = 0.01
z = 0
t = 0
dt = 1# / 100#
a(1) = 0#
a(2) = 2# / 27#
a(3) = 1# / 9#
a(4) = 1# / 6#
a(5) = 5# / 12#
a(6) = 0.5
a(7) = 5# / 6#
a(8) = 1# / 6#
a(9) = 2# / 3#
a(10) = 1# / 3#
a(11) = 1#
a(12) = 0#
a(13) = 1
For j = 1 To 13
For k = 1 To 13
b(j, k) = 0
Next k
kx(j) = 0#
ky(j) = 0#
kz(j) = 0#
Next j
b(2, 1) = 2# / 27#
b(3, 1) = 1# / 36
b(3, 2) = 1# / 12#
b(4, 1) = 1# / 24#
b(4, 3) = 1# / 8#
b(5, 1) = 5# / 12#
b(5, 3) = -25# / 16#
b(5, 4) = 25# / 16#
b(6, 1) = 1# / 20#
b(6, 4) = 1# / 4#
b(6, 5) = 1# / 5#
b(7, 1) = -25# / 108#
b(7, 4) = 125# / 108#
b(7, 5) = -65# / 27#
b(7, 6) = 125# / 54#
b(8, 1) = 31# / 300#
b(8, 5) = 61# / 225#
b(8, 6) = -2# / 9#
b(8, 7) = 13# / 900#
b(9, 1) = 2#
b(9, 4) = -53# / 6#
b(9, 5) = 704# / 45#
b(9, 6) = -107# / 9#
b(9, 7) = 67# / 90#
b(9, 8) = 3#
b(10, 1) = -91# / 108#
b(10, 4) = 23# / 108#
b(10, 5) = -976# / 135#
b(10, 6) = 311# / 54#
b(10, 7) = -19# / 60#
b(10, 8) = 17# / 6#
b(10, 9) = -1# / 12#
b(11, 1) = 2383# / 4100#
b(11, 4) = -341# / 164#
b(11, 5) = 4496# / 1025#
b(11, 6) = -301# / 82#
b(11, 7) = 2133# / 4100#
b(11, 8) = 45# / 82#
b(11, 9) = 45# / 164#
b(11, 10) = 18# / 41#
b(12, 1) = 3# / 205#
b(12, 6) = -6# / 41#
b(12, 7) = -3# / 205#
b(12, 8) = -3# / 41#
b(12, 9) = 3# / 41#
b(12, 10) = 6# / 41#
b(13, 1) = -1777# / 4100#
b(13, 4) = -341# / 164#
b(13, 5) = 4496# / 1025#
b(13, 6) = -289# / 82#
b(13, 7) = 2193# / 4100#
b(13, 8) = 51# / 82#
b(13, 9) = 33# / 164#
b(13, 10) = 12# / 41#
b(13, 12) = 1#
c(6) = 34# / 105#
c(7) = 9# / 35#
c(8) = 9# / 35#
c(9) = 9# / 280#
c(10) = 9# / 280#
c(12) = 41# / 840#
c(13) = 41# / 840#
For i = 0 To 10000
Worksheets("Sheet1").Cells(i + 1, 3) = t
Worksheets("Sheet1").Cells(i + 1, 4) = x
Worksheets("Sheet1").Cells(i + 1, 5) = y
Worksheets("Sheet1").Cells(i + 1, 6) = z
kx(1) = f1(t + dt, x, y, z) * dt
ky(1) = f2(t + dt, x, y, z) * dt
kz(1) = f3(t + dt, x, y, z) * dt
For j = 2 To 13
xd = x
yd = y
zd = z
For k = 1 To j - 1
xd = xd + b(j, k) * kx(k)
yd = yd + b(j, k) * ky(k)
zd = zd + b(j, k) * kz(k)
Next k
kx(j) = f1(t + a(j) * dt, xd, yd, zd) * dt
ky(j) = f2(t + a(j) * dt, xd, yd, zd) * dt
kz(j) = f3(t + a(j) * dt, xd, yd, zd) * dt
Next j
For j = 1 To 13
x = x + c(j) * kx(j)
y = y + c(j) * ky(j)
z = z + c(j) * kz(j)
Next j
t = t + dt
Next i
Application.ScreenUpdating = True
End Sub
Function f1(t As Double, x As Double, y As Double, z As Double) As Double
f1 = -s * (x - y)
End Function
Function f2(t As Double, x As Double, y As Double, z As Double) As Double
f2 = -y - x * z + r * x
End Function
Function f3(t As Double, x As Double, y As Double, z As Double) As Double
f3 = x * y - bb * z
End Function
« Fallen Physicist, Rising Engineer | トップページ | クルトガ(ハイグレード)を買った。& クルトガの特許 »
「学問・資格」カテゴリの記事
- 高周波(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)
« Fallen Physicist, Rising Engineer | トップページ | クルトガ(ハイグレード)を買った。& クルトガの特許 »
コメント