Loading [MathJax]/jax/element/mml/optable/GeneralPunctuation.js

[matplotlib animation] 18.カッシーニの卵形線

matplotlib Animation

はじめに

matplotlib FuncAnimationでカッシーニの卵形線のアニメーションを作成する。

コード

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from IPython.display import HTML
fig, ax = plt.subplots()
def update(num):
if len(tamagos) > 0:
tamagos.pop().remove()
tamagoss.pop().remove()
r_p_p = np.sqrt(a**2 * np.cos(2*t) + np.sqrt(a**4 * np.cos(2*t)**2 +k[num]**4 - a**4))
r_p_m = np.sqrt(a**2 * np.cos(2*t) - np.sqrt(a**4 * np.cos(2*t)**2 +k[num]**4 - a**4))
x_p_p = r_p_p * c
y_p_p = r_p_p * s
x_p_m = r_p_m * c
y_p_m = r_p_m * s
tamago_p, = plt.plot(x_p_p, y_p_p, 'go', lw=2)
tamago_m, = plt.plot(x_p_m, y_p_m, 'mo', lw=2)
tamagos.append(tamago_p)
tamagoss.append(tamago_m)
ax.set_title("cassini, parameter $a$= "+str(a)+", $k$= "+str(k[num])[:4]+"")
ax.grid()
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_aspect('equal')
#plot data
t = np.linspace(0, 2*np.pi, 500)
a=1
k= np.hstack([np.linspace(0.1, .8, 50),np.linspace(0.8, 1.2, 100),np.linspace(1.2, 2, 50)])
c, s = np.cos(t), np.sin(t)
tamagos=[]
tamagoss=[]
ani = animation.FuncAnimation(fig, update, 200, interval=100)
HTML(ani.to_html5_video())
#dpi=100
#ani.save('tamago.mp4', writer="ffmpeg",dpi=dpi)

解説

モジュールのインポート

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from IPython.display import HTML

グリッドと軸ラベル、アスペクト比の設定

ax.grid()
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_aspect('equal')

卵形線データの生成

t = np.linspace(0, 2*np.pi, 500)
a=1
k= np.hstack([np.linspace(0.1, .8, 50),np.linspace(0.8, 1.2, 100),np.linspace(1.2, 2, 50)])
c, s = np.cos(t), np.sin(t)
r_p_p = np.sqrt(a**2 * np.cos(2*t) + np.sqrt(a**4 * np.cos(2*t)**2 +k[num]**4 - a**4))
r_p_m = np.sqrt(a**2 * np.cos(2*t) - np.sqrt(a**4 * np.cos(2*t)**2 +k[num]**4 - a**4))
x_p_p = r_p_p * c
y_p_p = r_p_p * s
x_p_m = r_p_m * c
y_p_m = r_p_m * s

卵形線は↓の式のように定義される。

(x^2+y^2)^2 -2a^2(x^2–y^2) = k^4 – a^4

今回は、媒介変数(極座標)で卵形線はデータを生成する。
θ(t)を0〜2π,パラメータaをa=1として、
パラメータkを順次変化させたアニメーションを作成する。
cとsをそれぞれcos(t)とsin(t)としている。

卵形線を媒介変数表示すると↓のようになる。

r^2 = a^2cos2\theta \pm \sqrt{a^4cos^2 2\theta +k^4 – a^4}

rに±があるので√が正なのをr_p_pとし、√が負なのをr_p_mとしている。

それぞれのrでx,yのデータを得ている。

アニメーションを表示するための設定

tamagos=[]
tamagoss=[]

アニメーションで表示したいグラフを入れる空のリスト、tamagos,tamagossをつくる。

アニメーションの設定

def update(num):
if len(tamagos) > 0:
tamagos.pop().remove()
tamagoss.pop().remove()
r_p_p = np.sqrt(a**2 * np.cos(2*t) + np.sqrt(a**4 * np.cos(2*t)**2 +k[num]**4 - a**4))
r_p_m = np.sqrt(a**2 * np.cos(2*t) - np.sqrt(a**4 * np.cos(2*t)**2 +k[num]**4 - a**4))
x_p_p = r_p_p * c
y_p_p = r_p_p * s
x_p_m = r_p_m * c
y_p_m = r_p_m * s
tamago_p, = plt.plot(x_p_p, y_p_p, 'go', lw=2)
tamago_m, = plt.plot(x_p_m, y_p_m, 'mo', lw=2)
tamagos.append(tamago_p)
tamagoss.append(tamago_m)
ax.set_title("cassini, parameter $a$= "+str(a)+", $k$= "+str(k[num])[:4]+"")

アニメーションでは古いプロットを消す必要があるので、
グラフのリスト(tamagos,tamagoss)にグラフが入っている場合は、pop().remove()で消す。

r_p_p,r_p_mでkの値が異なるデータを順次作成し、それぞれのrでx,yデータを生成してプロットする。goで緑丸での表示となる。moでマゼンタ丸での表示となる。

作成したプロットを空のリスト(tamagos,tamagoss)に入れることでアニメーションを表示できる。

ax.set_title(“cassini, parameter $a$= “+str(a)+”, $k$= “+str(k[num])[:4]+””)でタイトルにアニメーションで変化するkの値をリアルタイムに表示する。[:4]とすることで小数点以下が大量に表示されるのを防いでいる。

アニメーションの表示

ani = animation.FuncAnimation(fig, update, 200, interval=100)
HTML(ani.to_html5_video())

FuncAnimationでanimationの表示。updateを200step実行してアニメーションとする。intervalは100 ms なので20秒のアニメーションとなる。
HTML(ani.to_html5_video())でjupyter notebook上にアニメーションを表示できる。

アニメーションの保存

#dpi=100
#ani.save('tamago.mp4', writer="ffmpeg",dpi=dpi)

ani.save(‘filename.mp4’, writer=”ffmpeg”,dpi=dpi)でアニメーションを保存することができる。dpiを設定することで指定の解像度での保存が可能。

コードをダウンロード(.pyファイル)

コードをダウンロード(.ipynbファイル)

参考

Cassini oval - Wikipedia

コメント