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