Pythonで（比較的）簡単に計算できる非線形微分方程式

はじめに

不規則な状態を実現するために、非線形微分方程式を解きたいと思ったので、Pythonで簡単に実装できるいくつかの非線形微分方程式をまとめておきます。

環境

• Windows 10 home

• Python(3.7.6)

• Numpy(1.19.4)

• Scipy(1.5.2)

• matplotlib(3.3.2)

• JupyterLab(2.2.8)

Lorenz63

\begin{align}
&\frac{dx}{dt} = \sigma(x - y) \\
&\frac{dy}{dt} = x(\rho - z) -y \\
&\frac{dz}{dt} = xy - \beta z
\end{align}

Lorenzが1963年に発表したカオス的なふるまいをする非線形常微分方程式です。パラメータは$\sigma$ , $\rho$, $\beta$です。Python Scipyのsolve_ivpを利用して、ルンゲクッタ法で計算します。

import numpy as np
from scipy.integrate import solve_ivp
#Lorenz63
def lorenz63(t, xyz, sigma=10, r=28, b=8/3.0):
    x = xyz[0]
    y = xyz[1]
    z = xyz[2]
    dxdt = sigma * (-x + y)
    dydt = -x * z + r * x - y
    dzdt = x * y - b * z
    return dxdt, dydt, dzdt

#数値計算
t = np.arange(0, 40, 0.01)
sol_l63 = solve_ivp(lorenz63, t_span=[t[0],t[-1]], y0=[1., 1., 1.], t_eval=t)

#計算結果
x = sol_l63.y[0]    
y = sol_l63.y[1]
z = sol_l63.y[2]

#計算結果可視化
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from IPython.display import HTML
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure()
ax = plt.axes(projection='3d')
ax.set_xlabel('$x$', fontsize=16)
ax.set_ylabel('$y$', fontsize=16)
ax.set_zlabel('$z$', fontsize=16)
ims = []
for i in range(x.size):
    if (i % 50 == 0) or (i == 0):
        im = ax.plot(x[:i+1], y[:i+1], z[:i+1], c='g')
        ims.append(im)
ani_l63 = animation.ArtistAnimation(fig, ims, interval=100)
HTML(ani_l63.to_html5_video())

Lorenz96

\begin{align}
&\mathrm{For} \ i=1,..., N:\\
&\frac{dx_i}{dt} = (x_{i+1} - x_{i-2})x_{i-1} - x_i + F
\end{align}

Lorenzが1996年に発表した非線形常微分方程式です。データ同化の検証として利用されているのを見ます。N個の変数を利用し、$x_{-1} = x_{N-1}$, $x_0=x_N$, $x_{N+1}=x_1$という関係を持っています。Fがパラメータで、F=8が一般的に使われるようです。solve_ivpを利用してルンゲクッタ法で計算します。

import numpy as np
from scipy.integrate import solve_ivp
#Lorenz96
def lorenz96(t, x, N=5, F=8):
    d = np.empty(N)
    for i in range(N):
        d[i] = (x[(i+1) % N] -x[i-2]) * x[i-1] - x[i] + F
    return d

#数値計算
N=5
F=8
x0 = F * np.ones(N)
x0[0] += 0.01
t = np.arange(0, 40., 0.01)
sol_l96 = solve_ivp(lorenz96, t_span=[t[0], t[-1]], y0=x0, t_eval=t, args=(N, F))

#計算結果
x = sol_l96.y[0]
y = sol_l96.y[1]
z = sol_l96.y[2]

#計算結果可視化
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from IPython.display import HTML
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure()
ax = plt.axes(projection='3d')
ax.set_xlabel('$x$', fontsize=16)
ax.set_ylabel('$y$', fontsize=16)
ax.set_zlabel('$z$', fontsize=16)
ims = []
for i in range(x.size):
    if (i % 50 == 0) or (i == 0):
        im = ax.plot(x[:i+1], y[:i+1], z[:i+1], c='g')
        ims.append(im)
ani_l96 = animation.ArtistAnimation(fig, ims, interval=100)
HTML(ani_l96.to_html5_video())

KdV方程式

\frac{\partial u}{\partial t} + \alpha u \frac{\partial u}{\partial x} + \beta \frac{\partial^3 u}{\partial x^3} = 0

ソリトンを求める非線形方程式の一つです。$\alpha$と$\beta$がパラメータです。偏微分方程式のため、空間方向と時間方向の解法を考える必要があります。空間方向はフーリエ・スペクトル法を、時間方向はsolve_ivpで陰解法ルンゲクッタを用いて計算します。フーリエ・スペクトル法は、フーリエ級数で$u$を展開し計算する方法です。フーリエ級数を利用しているため、周期境界条件になります。フーリエ・スペクトル法は、波数空間において、変数に$ik$(i:虚数単位, $k$:波数)かけると微分を表現できるので、フーリエ変換のライブラリがあれば、空間微分を簡単に書くことができます。実数のフーリエ変換のため、ライブラリはnumpy.fft.rfft, irfftを利用しています。

import numpy as np
from numpy.fft import rfft, irfft, fftfreq
from scipy.integrate import solve_ivp

#KdV方程式
def kdv(t, u, k, a, b):#k:波数, a,b:パラメータ
    v = rfft(u)#フーリエ変換
    ux = irfft(1.0j * k * v)#1階微分
    conv = u * ux 
    disp = irfft((1.0j * k)**3 * v)
    return -a * conv - b * disp

#条件設定
L = 2.0 #x軸長さ
T = 10 #時間
N = 512 #格子数
m = N//2 + 1 #波数の大きさ
dt = 1.0e-2 #タイムステップ
x = np.linspace(0, L, N) #x軸空間格子点
t = np.arange(0, T, dt) #時間
k = fftfreq(m) * m * 2 * np.pi / L #波数

#パラメータ
a = 1.0
b = 0.022**2
#初期値
u0 = np.cos(np.pi * x)

#数値計算
sol_kdv = solve_ivp(kdv, [0, T], u0, method='Radau', t_eval=t, args=(k, a, b))

#計算結果
u = sol_kdv.y

#可視化
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from IPython.display import HTML
fig = plt.figure()
ims = []
for uu in u.T[::10]:
    im = plt.plot(x, uu, c='g')
    plt.xlabel("$x$")
    plt.ylabel("$u$")
    ims.append(im)
ani_kdv = animation.ArtistAnimation(fig, ims, interval=50)
HTML(ani_kdv.to_html5_video())

1次元Burgers方程式

\frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} = \nu \frac{\partial^2 u}{\partial x^2}

Navier-Stokes方程式の圧力項を無視した移流拡散方程式です。数値流体計算の検証や圧縮流体で見ることがあります。$\nu$はパラメータで、流体の場合は動粘性係数を表します。これも偏微分方程式であるため、空間方向と時間方向の解法を考える必要があります。空間方向にチェビシェフ・スペクトル法を、時間方向はsolve_ivpで陰解法ルンゲクッタを用いて計算します。チェビシェフ・スペクトル法は$u$をチェビシェフ多項式で展開する方法です。チェビシェフ多項式は区間$[-1, 1]$ において、直交性があります。フーリエ・スペクトル法では周期境界条件でしたが、チェビシェフ・スペクトル法では、ディリクレ境界条件を与えることができ、この計算例の境界条件では、$u(-1)=0, u(1)=0$を与えています。Python Numpyにチェビシェフ多項式のライブラリがないので、微分演算子を作成していますが、空間微分を、微分演算子と変数の積をとるというシンプルな形で表現することができます。（参考:Spectral Method in Matlab）

import numpy as np
from scipy.integrate import solve_ivp

#演算子と格子点を出力する
def cheb(N):
    '''Chebyshev polynomial differentiation matrix.
       Ref.: Trefethen's 'Spectral Methods in MATLAB' book.
    '''
    x      = np.cos(np.pi*np.arange(0,N+1)/N)
    if N%2 == 0:
        x[N//2] = 0.0 # only when N is even!
    c      = np.ones(N+1); c[0] = 2.0; c[N] = 2.0
    c      = c * (-1.0)**np.arange(0,N+1)
    c      = c.reshape(N+1,1)
    X      = np.tile(x.reshape(N+1,1), (1,N+1))
    dX     = X - X.T
    D      = np.dot(c, 1.0/c.T) / (dX+np.eye(N+1))
    D      = D - np.diag( D.sum(axis=1) )
    return D,x

#Burgers方程式
def burgers(t, u, nu, D, D2):
    conv = u * D @ u #移流項: D:1階微分
    vis = nu * D2 @ u #粘性項: D2:2階微分
    bur = -conv + vis
    bur[0] = 0 #境界条件
    bur[-1] = 0 #境界条件
    return bur

N = 200 #格子点数
T = 1#最大時間
t = np.linspace(0, T, 50) #時間
D, x = cheb(N) #D:1階微分演算子、x:x座標
D2 = D @ D #2階微分演算子

#パラメータ
nu = 1.5e-2 #動粘性係数

#初期値
u0 = np.sin(np.pi * x * 2)

#数値計算
sol_bur = solve_ivp(burgers, [0, T], u0, method='Radau', t_eval=t, args=(nu,D,D2))

#計算結果
u = sol_bur.y

#可視化
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from IPython.display import HTML
fig = plt.figure()
ims = []
for uu in u.T:
    im = plt.plot(x, uu, c='g')
    plt.xlabel("$x$")
    plt.ylabel("$u$")
    ims.append(im)
ani_bur = animation.ArtistAnimation(fig, ims, interval=50)
HTML(ani_bur.to_html5_video())

参考

Lorenz, Edward, "Determinitic nonperiodic flow", Journal of the Atmospheric Science, 1963.

Lorenz, Edward, ”Predictability - A problem partly solved”、Seminar on Predictability, 1996.

Lloyd N. Trefethen , "Spectral Method in Matlab", Society for Industrial and Applied Mathematics, 2001.

石岡圭一, "スペクトル法による数値計算入門", 東京大学出版会, 2004.