PythonでのLorenz回転動画作成

Pythonを利用してLorenz回転動画を作成してみました。matplotlibで各フレームの画像を作成し、ffmpegで1分弱のMP4フォーマットに変換しています。ffmpegについては、pipではなく個別のインストールが必要となります。

サンプルソース

import matplotlib.pyplot as plt
import numpy as np
import matplotlib.animation as animation
import datetime

#------------------------------------------------------------------------------
# 変数設定
#------------------------------------------------------------------------------
#ローレンツ初期値
s=10
r=28
b=8/3

# 時間の設定
dt = 0.005
num_steps = 20000

# 描画領域設定
fig  = plt.figure(figsize=(16,9), dpi=256)
ax = fig.add_subplot(projection='3d')

#------------------------------------------------------------------------------
# ローレンツ計算
#------------------------------------------------------------------------------
def lorenz(x, y, z):
    x_dot = s*(y - x)
    y_dot = r*x - y - x*z
    z_dot = x*y - b*z
    return x_dot, y_dot, z_dot

#------------------------------------------------------------------------------
# ローレンツ描画
#------------------------------------------------------------------------------
def lorenzView():

    # 空配列生成
    xs = np.empty(num_steps + 1)
    ys = np.empty(num_steps + 1)
    zs = np.empty(num_steps + 1)

    # 初期条件
    xs[0]=1
    ys[0]=1
    zs[0]=1

    # 点列の座標生成
    for i in range(num_steps):
        x_dot, y_dot, z_dot = lorenz(xs[i], ys[i], zs[i])
        xs[i + 1] = xs[i] + (x_dot * dt)
        ys[i + 1] = ys[i] + (y_dot * dt)
        zs[i + 1] = zs[i] + (z_dot * dt)

    # 座標軸設定
    ax.view_init(elev=20, azim=110)
    ax.set_xlabel("x axis", size=12)
    ax.set_ylabel("y axis", size=12)
    ax.set_zlabel("z axis", size=12)
    ax.set_title("Lorenz", size=12)

    # プロット設定
    ax.plot(xs, ys, zs, lw=0.3)
    ax.set_box_aspect((1,1,1))
    plt.subplots_adjust(left=0, right=1, bottom=0, top=1)

#------------------------------------------------------------------------------
# アニメーション描画
#------------------------------------------------------------------------------
def animate(i):
    ax.view_init(elev=30., azim=-60 +3.6*i)
    return fig,

#------------------------------------------------------------------------------
#メイン関数
#------------------------------------------------------------------------------
def main():

    #--------------------------------
    #開始メッセージ
    #--------------------------------
    print("描画開始")
    dt_now = datetime.datetime.now()
    print(dt_now)

    #--------------------------------
    #ローレンツ描画
    #--------------------------------
    lorenzView()

    #--------------------------------
    #回転アニメーション描画
    #--------------------------------
    ani = animation.FuncAnimation(fig, animate,frames=400, interval=150, blit=True)
    ani.save('rotate_lorenz.mp4', writer="ffmpeg",dpi=256)

    #--------------------------------
    #終了メッセージ
    #--------------------------------
    print("描画終了")
    dt_now = datetime.datetime.now()
    print(dt_now)


if __name__ == "__main__":
    main()

ソースの解析

Lorenzの座標を計算します。

#------------------------------------------------------------------------------
# ローレンツ計算
#------------------------------------------------------------------------------
def lorenz(x, y, z):
    x_dot = s*(y - x)
    y_dot = r*x - y - x*z
    z_dot = x*y - b*z
    return x_dot, y_dot, z_dot

Lorenzの描画を行います。

#------------------------------------------------------------------------------
# ローレンツ描画
#------------------------------------------------------------------------------
def lorenzView():

    # 空配列生成
    xs = np.empty(num_steps + 1)
    ys = np.empty(num_steps + 1)
    zs = np.empty(num_steps + 1)

    # 初期条件
    xs[0]=1
    ys[0]=1
    zs[0]=1

    # 点列の座標生成
    for i in range(num_steps):
        x_dot, y_dot, z_dot = lorenz(xs[i], ys[i], zs[i])
        xs[i + 1] = xs[i] + (x_dot * dt)
        ys[i + 1] = ys[i] + (y_dot * dt)
        zs[i + 1] = zs[i] + (z_dot * dt)

    # 座標軸設定
    ax.view_init(elev=20, azim=110)
    ax.set_xlabel("x axis", size=12)
    ax.set_ylabel("y axis", size=12)
    ax.set_zlabel("z axis", size=12)
    ax.set_title("Lorenz", size=12)

    # プロット設定
    ax.plot(xs, ys, zs, lw=0.3)
    ax.set_box_aspect((1,1,1))
    plt.subplots_adjust(left=0, right=1, bottom=0, top=1)

フレーム毎に回転位置を変更させます。

#------------------------------------------------------------------------------
# アニメーション描画
#------------------------------------------------------------------------------
def animate(i):
    ax.view_init(elev=30., azim=-60 +3.6*i)
    return fig,

Lorenzの描画から回転までのアニメーションを作成します。

#------------------------------------------------------------------------------
#メイン関数
#------------------------------------------------------------------------------
def main():

    #--------------------------------
    #開始メッセージ
    #--------------------------------
    print("描画開始")
    dt_now = datetime.datetime.now()
    print(dt_now)

    #--------------------------------
    #ローレンツ描画
    #--------------------------------
    lorenzView()

    #--------------------------------
    #回転アニメーション描画
    #--------------------------------
    ani = animation.FuncAnimation(fig, animate,frames=400, interval=150, blit=True)
    ani.save('rotate_lorenz.mp4', writer="ffmpeg",dpi=256)

    #--------------------------------
    #終了メッセージ
    #--------------------------------
    print("描画終了")
    dt_now = datetime.datetime.now()
    print(dt_now)

作成画像

Maximaでトーラスの結び目を描画
目次 Maximaでのソースコード下記設定にて、メッシュの粒度を設定します。記設定にて、軸座標表示を非表示… 続きを読む »
Maxima plot3d ファイル出力
目次 Maximaでのソースコード（平面）出力画像関連記事最近の投稿関連記事: Maximaで波模… 続きを読む »
Maximaで世界地図を描画
目次 Maximaでのソースコード（平面）出力画像 Maximaでのソースコード（地球儀）出力画像関連記… 続きを読む »
Maximaで波模様を描画
目次 Maximaでのソースコード「 grid 」にて、メッシュの粒度を設定します。「 box 」にて、軸… 続きを読む »
Maximaで巻貝モデルを描画
目次 Maximaでのソースコード「 grid 」にて、メッシュの粒度を設定します。「 box 」にて、軸… 続きを読む »