Sympyでオイラーラグランジュ方程式を解く

自分用のメモのために、記事を書いてみました。

#はじめに
運動方程式を導出するために、オイラーラグランジュ方程式を解くことがあると思います。しかしながら多変数のオイラーラグランジュ方程式を解くとき、計算量が多く面倒のなので無料で利用できるPythonの数式処理ライブラリであるSympyを使ってオイラーラグランジュ方程式を解いてみることにしました。

今回説明すること

一般的な運動である１質点とばねによる振動系のラグランジアンからオイラーラグランジュ方程式の解くプログラムを用いて説明をします。

コード

Euler–Lagrange.py

#Sympyをインポート
import sympy as sp
#JupyterNotebookを使っている場合数式を自然数表示
sp.init_printing()

#時間の変数を設定
t = sp.symbols('t')
#質点の位置の変数を設定
#（xはｔの関数に設定するためにオプションにcls=sp.Functionを追加）
x = sp.symbols('x ', cls=sp.Function)
#質量とばね定数の設定
m, k = sp.symbols('m k')

#ラグランジアンLを設定
L =sp.Function("L")
L = (m*(x(t).diff(t))**2)/2 -(k*x(t)**2)/2

#変数の座標と座標の時間微分に配列を作成
#多変数に拡張する場合ここに変数を追加
pos = [x(t)]
vel =[x(t).diff(t)]

#オイラーラグランジュ方程式解く関数を定義
def EulerLagrange(L,pos,vel):
    
    EQ_list =[]

    for i in range(len(pos)):
        EQ_list.append(sp.simplify(sp.Eq(L.diff(vel[i]).diff(t) - L.diff(pos[i]),0)))
    return EQ_list

#運動方程式
f =EulerLagrange(L,pos,vel)[0]
print(f)

上のように計算を行わせると、運動方程式$f$は、

k x{\left (t \right )} + m \frac{d^{2}}{d t^{2}} x{\left (t \right )} = 0

となり、１質点とばねによる振動系と同様の形になりました。

これで、自宅でも簡単に運動方程式が導出できます。

#終わりに
運動方程式の導出ができました。実際に欲しいの運動方程式を数値計算した結果が欲しかったのですが、この求めた運動方程式をPythonのライブラリ（Sympy、Numpy、SciPyなど）を用いて数値計算する方法がわかりませんせんでした。もしご存じの方がいらっしゃいましたら、教えていたただけると幸いです。