« 餃子の王将で7月度限定メニュー、酸辣湯麺(Bセット)を食す。そのままでも美味しいが、卓上のお酢をたっぷりかけて辣油もたらすとさらにおいしく。 | トップページ | 松屋で店舗限定のめちゃ大盛麻婆豆腐丼を食す。確かに餡の量がはんぱない。 »

2020年7月16日 (木)

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アニメが始まります。(最初が動きがないので面白くないがだんだん動く)

Kuramoto_sivashinsky1

 

ソースは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セット)を食す。そのままでも美味しいが、卓上のお酢をたっぷりかけて辣油もたらすとさらにおいしく。 | トップページ | 松屋で店舗限定のめちゃ大盛麻婆豆腐丼を食す。確かに餡の量がはんぱない。 »

学問・資格」カテゴリの記事

コメント

コメントを書く

(ウェブ上には掲載しません)

« 餃子の王将で7月度限定メニュー、酸辣湯麺(Bセット)を食す。そのままでも美味しいが、卓上のお酢をたっぷりかけて辣油もたらすとさらにおいしく。 | トップページ | 松屋で店舗限定のめちゃ大盛麻婆豆腐丼を食す。確かに餡の量がはんぱない。 »

最近の記事

2021年5月
            1
2 3 4 5 6 7 8
9 10 11 12 13 14 15
16 17 18 19 20 21 22
23 24 25 26 27 28 29
30 31          
フォト
無料ブログはココログ