Python+Scipyでルンゲクッタ8次のDOP853(Dormand&Prince)を使う(その14) 蔵本・シバシンスキー方程式(Kuramoto-Sivashinsky equation)を計算してGIFアニメに。
蔵本先生と言えば、最近は千葉先生の研究で蔵本モデルが有名ですが、私にはこっちのほうがなじみがある蔵本・シバシンスキー方程式(Kuramoto-Sivashinsky equation)。1次元の乱流がシンプルに見られる。これはKdV方程式ができればすぐにできる。
∂u/∂t +u*∂u/∂x = -∂^2u/∂x^2 -∂^4/∂x^4
空間微分をO(h^6)という高次の差分で近似、かつ時間微分をルンゲクッタ8次のスキームであるDormand-Prince法(DOP853)で計算。
結果はこんな感じ。クリックでGIFアニメが始まります。(最初が動きがないので面白くないがだんだん動く)
ソースはKdVと同じなので関数部分だけ。
def kuramoto(t, u): #odeintのときとt,xの並びが逆
f=np.zeros(Nx)
for i in range(Nx):
if i < Nx-1:
p1 = i + 1
else:
p1 = 0
if i < Nx - 2:
p2 = i + 2
elif i == Nx - 2:
p2 = 0
else:
p2 = 1
if i < Nx - 3:
p3 = i + 3
elif i == Nx - 3:
p3 = 0
elif i == Nx - 2:
p3 = 1
else:
p3 = 2
if i > 0:
m1 = i - 1
else:
m1 = Nx-1
if i > 1:
m2 = i - 2
elif i == 1:
m2 = Nx - 1
else:
m2 = Nx - 2
if i > 2:
m3 = i - 3
elif i == 2:
m3 = Nx - 1
elif i == 1:
m3 = Nx - 2
else:
m3 = Nx - 3
d1 = (u[p3] / 9 - u[p2] + 5* u[p1] - 5 * u[m1] + u[m2] - u[m3] / 9) / (20 * dx / 3)
d2 = (2 * u[p3] - 27 * u[p2] + 270 * u[p1] - 490 * u[i] + 270 * u[m1] - 27 * u[m2] + 2 * u[m3])/(180 * dx * dx)
d4 = (-u[p3] +12*u[p2]-39*u[p1]+56*u[i]-39*u[m1]+12*u[m2]-u[m3])/(6*dx*dx*dx*dx)
f[i] = -u[i] * d1 - d2 - d4
return f
« 餃子の王将で7月度限定メニュー、酸辣湯麺(Bセット)を食す。そのままでも美味しいが、卓上のお酢をたっぷりかけて辣油もたらすとさらにおいしく。 | トップページ | 松屋で店舗限定のめちゃ大盛麻婆豆腐丼を食す。確かに餡の量がはんぱない。 »
「学問・資格」カテゴリの記事
- 高周波(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)
« 餃子の王将で7月度限定メニュー、酸辣湯麺(Bセット)を食す。そのままでも美味しいが、卓上のお酢をたっぷりかけて辣油もたらすとさらにおいしく。 | トップページ | 松屋で店舗限定のめちゃ大盛麻婆豆腐丼を食す。確かに餡の量がはんぱない。 »
コメント