数学・物理学の具体的な計算にPythonを使って、数学・物理学もPythonも同時に学んでしまいましょう。前回の記事
では、Pythonを使って振動を表す運動方程式（2階常微分方程式）の解を数値的に求めました。今回の記事では、その数値解によって記述される質点の運動の様子をアニメーションによって可視化してみたいと思います。
本記事は以下の記事の続編的な位置付けになります。あわせて読んでみてください：
pianofisica.hatenablog.com
pianofisica.hatenablog.com
単振動の運動方程式とその数値解
単振動（１次元調和振動子）の運動方程式は以下で与えられます：
この方程式をどのようにして数値的に解くかについては、上記の過去記事の解説を参照してください。ここでは必要なコードだけ載せておきます：
N = 3000 # ステップ数 h = 0.01 # 差分間隔 X = [0] # 初期時刻 Y = [1] # 初期位置 P = [0] # 初期速度 # 差分方程式 def x(k): return h+X[k-1] def y(k): return Y[k-1]+h*(P[k-1])+h**2/2*(-Y[k-1]) def p(k): return P[k-1]+h*(-Y[k-1])+h**2/2*(-P[k-1]) for k in range(N): X.append( x(k+1) ) Y.append( y(k+1) ) P.append( p(k+1) ) import matplotlib.pyplot as plt # 図の描画 plt.plot(X, Y) # 図の表示 plt.show()
以上のプログラムから、質点の運動の様子として得られるのが次の図です：
ここで横軸は時刻、縦軸が質点の位置です。
アニメーションの作成
上で得られた質点の位置の時間的変化をアニメーションによって可視化してみましょう。アニメーションの作成にはmatplotlibライブラリの関数を利用します：
%matplotlib notebook # JupyterNotebookでアニメーションを表示させるために必要です import matplotlib.animation as animation fig = plt.figure() # グラフの設定 plt.axes().set_aspect('equal') plt.xlim(-1.5, 1.5) # x軸の表示範囲 plt.ylim(-0.5, 0.5) # y軸の表示範囲 # アニメーションのフレームの作成 frames = [] num_frames = int(N/10) for i in range(num_frames): x = Y[10*i] # アニメーションのためにプロットするデータ点は10個おき y = 0 frame = plt.scatter(x, y, color='blue') frames.append([frame]) # アニメーションの描画 ani = animation.ArtistAnimation(fig, frames, interval=50) # アニメーションの表示 plt.show() # GIFファイルとして保存 ani.save('harmonic_ani.gif', writer='pillow')
以上から次のようなアニメーションが出力されるとともに、ローカルに harmonic_ani.gif というファイル名のGIFファイルが作成されるはずです：
（追記）不用意にバージョンを更新しPythonの環境を変更したところ、Jupyter Notebookでアニメーションが出力できなくなってしまいました。一応の解決方法としては、一旦アニメーションのGIFを作成してそれを読み込ませる、というものです：
%matplotlib inline import matplotlib.animation as animation fig = plt.figure() # グラフの設定 plt.axes().set_aspect('equal') plt.xlim(-1.5, 1.5) # x軸の表示範囲 plt.ylim(-0.5, 0.5) # y軸の表示範囲 # アニメーションのフレームの作成 frames = [] num_frames = int(N/10) for i in range(num_frames): x = Y[10*i] # アニメーションのためにプロットするデータ点は10個おき y = 0 frame = plt.scatter(x, y, color='blue') frames.append([frame]) # アニメーションの描画 ani = animation.ArtistAnimation(fig, frames, interval=50) # アニメーションの表示 from IPython.display import Image f = r'harmonic_ani.gif' writergif = animation.PillowWriter(fps=50) ani.save(f,writer=writergif) plt.close() Image(open('harmonic_ani.gif','rb').read())
前回の記事（微分方程式の数値的解法）と今回の記事（アニメーションによる視覚化）で扱った内容によって、解析的な手法に頼ることなくさまざまな運動の様子について一定程度まで調べることができるようになりますね。
