<jupyter_start><jupyter_text>ソルバーを使わずに連立常微分方程式を解くここでは多次元の常微分方程式を扱う．  題材は有名なローレンツ方程式とする．     -->$$\begin{bmatrix}\dot{x}\\\dot{y}\\\dot{z}\\\end{bmatrix} = \begin{bmatrix}\ -px+py\\\ -xz+rx-y\\\ xy-bz\\\end{bmatrix}$$ただし  $$\begin{bmatrix}\ x_0\\\ y_0\\\ z_0\\\end{bmatrix} = \begin{bmatrix}\ 0\\\ 4\\\ 28\\\end{bmatrix},  \quad  p = 10,  \quad  r = 28,  \quad  b = \frac{8}{3}$$とする．<jupyter_code>import matplotlib.pyplot as plt
import numpy as np<jupyter_output><empty_output><jupyter_text>標準のリストで多次元配列の和を計算するのはかなり面倒です（実装してみるとわかる）．  ex. 次の２つのベクトル$a$と$b$の和を計算する  $$a = \begin{bmatrix}1\\2\\\end{bmatrix},~~~b = \begin{bmatrix}3\\4\\\end{bmatrix}$$<jupyter_code># 標準のリストを使う
a = [1, 2]
b = [3, 4]

ans = []
for i in range(2):
    ans.append(a[i] + b[i])

ans
# numpyで実装
a = np.array([1, 2])
b = np.array([3, 4])

a + b<jupyter_output><empty_output><jupyter_text>このようにnumpyを使うことで直感的な行列（ベクトル）演算が行なえます．  多次元の場合，ベクトル和を一気に計算できるのでndarrayを使うのが楽です．   微分方程式の定義<jupyter_code>def loorenz(X, p, r, b):
    """ローレンツ方程式
    * Xは状態変数ベクトル（ndarray）
    * p, r, bはパラメータ
    """
    
    x, y, z = X  # アンパック．右辺のndarryから１つずつ値を取り出す．
    
    dx = -p*x + p*y
    dy = -x*z + r*x -y
    dz = x*y - b*z
    
    return np.array([dx, dy, dz])<jupyter_output><empty_output><jupyter_text>ルンゲクッタ法で解いてみる．<jupyter_code>dt = 0.01

# パラメータ
p = 10
r = 28
b = 8/3

N_TOTAL = 5000  #ループ回数
x = []  # 解xを格納する空のリスト
t = []  # 時刻tnを格納するリスト

# 初期値代入
x.append(np.array([0, 4, 28]))
t.append(0)

for i in range(N_TOTAL-1):
    # 係数を計算
    k1 = loorenz(x[i], p, r, b)
    k2 = loorenz(x[i] + dt/2*k1, p, r, b)
    k3 = loorenz(x[i] + dt/2*k2, p, r, b)
    k4 = loorenz(x[i] + dt*k3, p, r, b)
    
    x_next = x[i] + dt/6*(k1 + 2*k2 + 2*k3 + k4)  # 解を更新
    t_next = t[i] + dt
    
    # 解を格納
    x.append(x_next)
    t.append(t_next)<jupyter_output><empty_output><jupyter_text>解`x`はndarrayのリストになっているので特定の行を取り出すのが大変です．<jupyter_code>x[:5]<jupyter_output><empty_output><jupyter_text>データにアクセスしやすくするため，単一のndarrayに変換する．<jupyter_code>x = np.stack(x)
x<jupyter_output><empty_output><jupyter_text>0列目にxの時刻歴，1列目にyの時刻歴，2列目にzの時刻歴が入っている二次元ndarray`x`が作れた．  解をプロットしてみる．<jupyter_code>plt.plot(t, x[:, 0], label="x")
plt.plot(t, x[:, 1], label="y")
plt.plot(t, x[:, 2], label="z")
plt.xlabel("t")
plt.legend()<jupyter_output><empty_output><jupyter_text>周期的な振動が起こっていることがわかる．  3次元プロットしてみる．<jupyter_code>fig = plt.figure(figsize=(8, 8))
ax = fig.add_subplot(projection='3d')
ax.plot(x[:, 0], x[:, 1], x[:, 2])
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_zlabel("z")<jupyter_output><empty_output><jupyter_text>有名なローレンツアトラクタの図が得られた．  jupyter以外で`plt.show()`すると別ウィンドウで3Dグラフが表示され，くるくる回転させたりズームさせたりできる．  ---- ndarrayに書き込む方法空リストに計算結果をappendしていくのではなく，はじめに結果を書き出すndarrayを準備してから順に書き込んでいくやり方です．  ~~pythonではリストにappendすると，前のリストを破棄して新たなリストを作成するため遅くなります．~~  ~~そのため予めndarray等で場所を確保してから書き込んでいくほうが高速です．~~<jupyter_code>x_nd = np.empty((N_TOTAL, 3))  # 解xを格納する二次元ndarray
t_nd = np.empty(N_TOTAL)  # 時刻tnを格納する一次元ndarray

# 初期値を書き込む
x_nd[0, :] = np.array([0, 4, 28])  # :はその次元の全てを選択することを表す．(0:end と等価)
t_nd[0] = 0

for i in range(N_TOTAL-1):
    k1 = loorenz(x_nd[i,:], p, r, b)
    k2 = loorenz(x_nd[i,:] + dt/2*k1, p, r, b)
    k3 = loorenz(x_nd[i,:] + dt/2*k2, p, r, b)
    k4 = loorenz(x_nd[i,:] + dt*k3, p, r, b)
    x_next = x_nd[i, :] + dt/6*(k1 + 2*k2 + 2*k3 + k4)  # 解を更新
    t_next = t[i] + dt
    # 解を格納
    x_nd[i+1, :] = x_next
    t_nd[i+1] = t_next
fig_nd = plt.figure(figsize=(8, 8))
ax_nd = fig_nd.add_subplot(projection='3d')
ax_nd.plot(x_nd[:, 0], x_nd[:, 1], x_nd[:, 2])
ax_nd.set_xlabel("x")
ax_nd.set_ylabel("y")
ax_nd.set_zlabel("z")<jupyter_output><empty_output><jupyter_text>~~ループ回数が1000回くらいでは違いは出ませんが，ループ回数が増えれば増えるほど計算時間の差は大きくなります．~~  jupyterのマジックコマンド`%timeit`を使って比較してみる．   速度比較<jupyter_code>def hutsu_no_list(N):
    x = []
    t = []
    x.append(np.array([0, 4, 28]))
    t.append(0)

    for i in range(N-1):
        k1 = loorenz(x[i], p, r, b)
        k2 = loorenz(x[i] + dt/2*k1, p, r, b)
        k3 = loorenz(x[i] + dt/2*k2, p, r, b)
        k4 = loorenz(x[i] + dt*k3, p, r, b)        
        x_next = x[i] + dt/6*(k1 + 2*k2 + 2*k3 + k4)
        t_next = t[i] + dt
        
        x.append(x_next)
        t.append(t_next)
    
    return t, np.stack(x)


def ndarray_tsukau(N):
    x_nd = np.empty((N, 3))
    t_nd = np.empty(N)
    x_nd[0, :] = np.array([0, 4, 28])
    t_nd[0] = 0
    

    for i in range(N-1):
        k1 = loorenz(x_nd[i,:], p, r, b)
        k2 = loorenz(x_nd[i,:] + dt/2*k1, p, r, b)
        k3 = loorenz(x_nd[i,:] + dt/2*k2, p, r, b)
        k4 = loorenz(x_nd[i,:] + dt*k3, p, r, b)
        x_nd[i+1, :] = x_nd[i, :] + dt/6*(k1 + 2*k2 + 2*k3 + k4)
        t_nd[i+1] = t_nd[i] + dt

    return t_nd, x_nd<jupyter_output><empty_output><jupyter_text>ループ回数を指定<jupyter_code>N = 100000
%%timeit
hutsu_no_list(N)
%%timeit
ndarray_tsukau(N)<jupyter_output>1.04 s ± 4.6 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)<jupyter_text>ndarrayを準備するほうが遅い？  numbaで高速化してみる．<jupyter_code>from numba import njit

@njit("void(float64[:], float64, float64, float64, float64[:])")
def loorenz_2(X, p, r, b, out):
    out[0] = -p*X[0] + p*X[1]
    out[1] = -X[0]*X[2] + r*X[0] - X[1]
    out[2] = X[0]*X[1] - b*X[2]

@njit("Tuple((float64[:], float64[:, :]))(int64)")
def ndarray_tsukau_2(N):
    x_nd = np.empty((N, 3))
    t_nd = np.empty(N)
    x_nd[0, :] = np.array([0, 4, 28])
    t_nd[0] = 0
    
    k1 = np.empty(3)
    k2 = np.empty(3)
    k3 = np.empty(3)
    k4 = np.empty(3)

    for i in range(N-1):
        loorenz_2(x_nd[i,:], p, r, b, k1)
        loorenz_2(x_nd[i,:] + dt/2*k1, p, r, b, k2)
        loorenz_2(x_nd[i,:] + dt/2*k2, p, r, b, k3)
        loorenz_2(x_nd[i,:] + dt*k3, p, r, b, k4)
        x_nd[i+1, :] = x_nd[i, :] + dt/6*(k1 + 2*k2 + 2*k3 + k4)
        t_nd[i+1] = t_nd[i] + dt

    return t_nd, x_nd
%%timeit
ndarray_tsukau_2(N)<jupyter_output>15 ms ± 10.2 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
{'fraction_non_alphanumeric': 0.20003390405153415, 'fraction_numerical': 0.043566706221393456, 'mean_word_length': 3.580745341614907, 'pattern_counts': {'":': 0, '<': 47, '<?xml version=': 0, '>': 48, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '12', 'dir': 'jupyter-structured-clean-dedup', 'id': '252408', 'license': 'MIT', 'n_tokens_mistral': 3617, 'n_tokens_neox': 3460, 'n_words': 485, 'path': 'ode/python_src/simple2.ipynb', 'repo_name': 'YoshimitsuMatsutaIe/abc_2022', 'size': '5899'}