力学

【Pythonで実装】N重振り子の運動方程式を解く。アニメーションあり

こんにちは(@t_kun_kamakiri)(^^)/

今回はN重振り子の運動方程式の問題を取り扱います。

本記事の内容

N重振り子の運動方程式を導く。

  • N重振り子の運動
  • ラグランジュ方程式から運動方程式を導く
  • N重振り子の数値計算

こちらの記事を読んで

振り子の運動を数値計算で再現したい方には参考になる記事です。

スポンサーリンク

N重振り子の運動方程式

本記事では一般的なN重振り子の問題を考えますが、慣れない方は2重振り子や3重振り子など自由度の少ないの運動方程式から考えてみても良いでしょう。
一般的には自由度が大きくなるとそれだけ問題が複雑になるのですが、逆に自由度が多いほうが帰って式がシンプルになる場合があります

なので、今回はあえてN個の質点が連なった振り子の問題を扱おうと思います。

  • 質点の質量:$m_{i}$
  • 質点の位置:$(x_{i},y_{i})$
  • 質点間の長さ:$l_{i}$

質点には重力加速度$m_{i}g$がはたらいています。

運動方程式と等価な方程式

力学の問題を解く際には、まず運動方程式である偏微分方程式を導出する必要があります。
運動方程式と等価となる方程式はいくつか方法かあって、問題に応じて適した方法を用いれば良いと思います。
代表的な運動方程式と等価な方程式を示します。

ニュートンの第二法則

\begin{align*}
ma=F
\end{align*}

直交する各成分ごとで運動方程式を立てて微分方程式を解く。※デカルト座標と極座標で運方程式の形が変わる。
※qの2回微分が出てくるので取り扱いがちょっと面倒

振り子の問題はデカルト座標で解くと変数が$2n$個になるため、極座標で解くと変数が$\theta$のみになり問題が簡単になります。運動方程式を導くまでにベクトルとしてい扱う必要がありどうしても導出が複雑になってしまいます。

オイラーラグランジュ方程式

\begin{align*}
\frac{d}{dt}\frac{\partial L}{\partial \dot{q}}=\frac{\partial L}{\partial q}
\end{align*}

一般化座標$q$と一般化速度$\dot{q}$を使って統一した一本の式で記述できるので、デカルト座標と極座標とでオイラーラグランジュ方程式の形は変わらない。
※qの2回微分が出てくるので取り扱いがちょっと面倒

ハミルトンの正準方程式

\begin{align*}
\dot{q}&=\frac{\partial H}{\partial p}\\
p&=-\frac{\partial H}{\partial q}
\end{align*}

一般化座標$q$と一般化運動量$p$を使って統一した一本の式で記述できます。
※$q$の2回微分が出てこない代わりに、変数の数が$p,q$の2つになり式の数が増えますが時間微分項が無い分問題は解きやすくなります。

このどれを使っても良いのですが、今回は簡単に運動方程式を導くことができるオイラー・ラグランジュ方程式を使った運動方程式の導出方法を紹介します

運動方程式を導く

運動方程式を導くを導くまでのフローを以下に記します。

\begin{align*}
x_{i}&=\sum_{j=1}l_{j}\sin\theta_{j}\tag{1}\\
y_{i}&=\sum_{j=1}-l_{j}\cos\theta_{j}\tag{2}
\end{align*}

まずは運動方程式を書き下します。

\begin{align*}
T&=\sum_{i=1}^{n}\frac{1}{2}m_{i}(\dot{x}_{i}^2+\dot{y}_{i}^2)\\
&=\sum_{j=1}^{n}\sum_{k=1}^{i}l_{j}l_{k}\dot{\theta}_{j}\dot{\theta}_{k}\cos(\theta_{j}-\theta_{k})\tag{3}
\end{align*}

次に位置エネルギーを書きます。

\begin{align*}
U&=\sum_{i=1}^{n}mgy_{i}\\
&=\sum_{i=1}^{n}mg\sum_{j=1}^{i}-l_{i}\cos\theta_{i}\tag{4}
\end{align*}

ここでの変数は$x_{i},y_{i}$ではなく、$l_{i},\theta_{i}$が、$l_{i}$は固定値なので実際取り扱う変数は$\theta_{i}$のみです。
そして、$N$個の質点があるので自由度は$N$ということになりますね。

ラグランジアンを求める

運動方程式を導くには、ラグランジアン$L=T-U$を求めてオイラー・ラグランジュ方程式に代入すれば良いので、まずはラグランジアンを求めます。

ラグランジュ方程式の導出の仕方は↓こちらの記事に書いていますのでご参考ください。

\begin{align*}
L&=T-U\\
&=\sum_{j=1}^{n}\sum_{k=1}^{i}l_{j}l_{k}\dot{\theta}_{j}\dot{\theta}_{k}\cos(\theta_{j}-\theta_{k})+\sum_{i=1}^{n}mg\sum_{j=1}^{i}l_{i}\cos\theta_{i}\tag{5}
\end{align*}

ここでの一般化座標は$\theta$だし、一般化速度は$\dot{\theta}$であることに注意して、次はオイラー・ラグランジュ方程式から運動方程式を導く計算を進めていきます。

オイラー・ラグランジュ方程式から運動方程式を導く

オイラー・ラグランジュ方程式は以下で与えられます。

\begin{align*}
\frac{d}{dt}\frac{\partial L}{\partial \dot{q}}=\frac{\partial L}{\partial q}\tag{6}
\end{align*}

(6)に(5)を代入します。

\begin{align*}
l_{i}\sum_{j=1}^{n}\bigg(\sum_{k=max(i,j)}^{n}m_{k}\bigg)l_{j}\left \{ \ddot{\theta_{j}}\cos(\theta_{i}-\theta_{j}) -\dot{\theta}_{j}\big(\theta_{i}-\theta_{j}\big)\sin(\theta_{i}-\theta_{j})\right \}\\
=-l_{i}\left \{\sum_{j=1}^{n}\bigg(\sum_{k=max(i,j)}^{n}m_{k}\bigg)\big(l_{j}\dot{\theta}_{i}\dot{\theta}_{j}\sin(\theta_{i}-\theta_{j})\big) +\big(\sum_{k=i}^{n}g\sin\theta_{i}\big) \right \}\tag{7}
\end{align*}

今回の運動方程式(偏微分方程式)における変数は$\theta$なので、$\theta$に関する2階微分方程式であることに注意します。

ここで、自由度が$N$個あるので方程式を整理すると$N$行$N$列の行列Aと$B$を使って整理することができます。

行列計算を解くだけということになります。

\begin{align*}
\ddot{\theta}=-A^{-1}BE\tag{8}
\end{align*}

行列$B$については、対角成分と非対角成分とで式の形が異なるので注意が必要です。

$N$重振り子の場合の運動方程式が(8)の形で書け、行列計算を行えば時刻$t$での$\theta(t)_{1},\theta(t)_{2},\cdots ,\theta(t)_{n}$がわかるということですね。

計算の途中で、逆行列$A^{-1}$を計算する必要がありますが、プログラミング言語Pythonの数値計算用のライブラリであるnumpyを使えば簡単に計算ができます。

PythonでN重振り子の運動の数値計算

Pythonを使って(8)の数値計算をしたいと思います。

今回は、以下のようにExcelに振り子のパラメータになる重量と質点間距離を書き込んでおき、PythonからExcelデータを読み込むプログラムにしたいと思います。
その方が、もしPythonを知らない人でもExcel内のパラメータを変えることで計算ができるようになるからです。

解くべき方程式である(8)$\ddot{\theta}=-A^{-1}BE$は$\theta$に関する2階微分方程式であるため、時間の積分を数値的に行う必要があります。

時間の数値積分は精度の良い4次のツンゲクッタを使って解くことにします。

jupyter labを使ってプログラムを作成します。

pandasからExcelを読み込む部分で XLRDError: Excel xlsx file; not supported というエラーになってしまいます。理由は、pd.read_excel() のデフォルトエンジンが xlrd のためでです。対策は pip install openpyxl で openpyxl をインストールしてpd.read_excel('N_param0.xlsx',header=1, engine='openpyxl')と書けば良いです。

まずはExcelを用意します。
質点の数は何個でも良いのですが、今回は5個にしておきます。
質点の数を増やしたい場合はG列以降にも値を書き込んでいけば自動的に質点の数が増えて振り子の計算してくれます。

必要なライブラリをインポートします。

次にExcelのデータをpandasで読み込みます。

ちゃんと読み込めているか確認しましょう。

値が正しく読み取れているか確認しておきます。

良さそうだったら、振り子の初期状態を確認しましょう。

微分方程式を解くプログラムを作成します。

後は、時間に対してどんどん積分を繰り返します。

ここまでで$\theta(t)_{1},\theta(t)_{2},\cdots,\theta(t)_{5}$が求まりましたが、振り子のアニメーションを作成するために、$x(t),y(t)$に変換する必要があります。

アニメーションにする前にどういった挙動をしているのかをグラフにしてみます。

正しそうな挙動をしていますね。

では、アニメーションを作成します。

 

うまくいったのではないでしょうか(^^♪

まとめ

N重振り子の運動方程式をオイラー・ラグランジュ方程式を使って導出し、Pythonで実装してアニメーションにしました。

 

【プロフィール】

カマキリ
(^^)

大学の専攻は物性理論で、Fortranを使って数値計算をしていました。
CAEを用いた流体解析は興味がありOpenFOAMを使って勉強しています。

プロフィール記事はこちら

 

大学学部レベルの物理の解説をします 大学初学者で物理にお困りの方にわかりやすく解説します。

このブログでは主に大学以上の物理を勉強して記事にわかりやすくまとめていきます。

  • ・解析力学
  • ・流体力学
  • ・熱力学
  • ・量子統計
  • ・CAE解析(流体解析)
  • note
    noteで内容は主に「プログラミング言語」の勉強の進捗を日々書いています。また、「現在勉強中の内容」「日々思ったこと」も日記代わりに書き記しています。
  • youtube
    youtubeではオープンソースの流体解析、構造解析、1DCAEの操作方法などを動画にしています。
    (音声はありません_(._.)_)
  • Qiita
    Qiitaではプログラミング言語の基本的な内容をまとめています。