座標の取り方は下図のように取ります。棒の伸び縮みは無いものとします。

どういう解き方でもいいですが、ここでは

1. デカルト座標\(L(x,y)\)でラグランジアンを記述
2. デカルト座標から座標変換し、\((r,\theta)\)でラグランジアンを記述
3. 新たな座標系で運動方程式を導く

という順で解いていきます。

[adsense1]

1, デカルト座標でのラグランジアンLは(運動エネルギーK)-(位置エネルギーU)と書けるため、
\(
L(x_1,\dot{x}_1,y_1,\dot{y}_1,x_2,\dot{x}_2,y_2,\dot{y}_2)\\
\displaystyle =\frac{1}{2}m_1(\dot{x}_1^2+\dot{y}_1^2)+\frac{1}{2}m_1(\dot{x}_1^2+\dot{y}_2^2)-(-m_1gy_1-m_2gy_2)
\)
と書けます。

2, デカルト座標から座標変換
式を簡単にするために座標変換を行います。新しい座標\((r_1,\theta_1,r_2,\theta_2)\)とデカルト座標\((x_1,y_1,x_2,y_2)\)の関係式は
\(
\begin{align}
x_1&=r_1\sin{\theta_1}\\
y_1&=-r_1\cos{\theta_1}\\
x_2&=r_1\sin{\theta_1}+r_2\sin{\theta_2}\\
y_2&=-r_1\cos{\theta_1}-r_2\cos{\theta_2}
\end{align}
\)
という関係があります。各々を時間で微分すれば、
\(
\begin{align}
\dot{x}_1&=\dot{r}_1\sin{\theta_1}+r_1\dot{\theta}_1\cos{\theta_1}\\
\dot{y}_1&=-\dot{r}_1\cos{\theta_1}+r_1\dot{\theta}_1\sin{\theta_1}\\
\dot{x}_2&=\dot{r}_1\sin{\theta_1}+r_1\dot{\theta}_1\cos{\theta_1}+\dot{r}_2\sin{\theta_2}+r_2\dot{\theta}_2\cos{\theta_2}\\
\dot{y}_2&=-\dot{r}_1\cos{\theta_1}+r_1\dot{\theta}_1\sin{\theta_1}-\dot{r}_2\cos{\theta_2}+r_2\dot{\theta}_2\sin{\theta_2}
\end{align}
\)
これらをラグランジアン\(L(x_1,\dot{x}_1,y_1,\dot{y}_1,x_2,\dot{x}_2,y_2,\dot{y}_2)\)に代入します。すると、新たな座標系でのラグランジアン\(L(r_1,\dot{r}_1,\theta_1,\dot{\theta}_1,r_2,\dot{r}_2,\theta_2,\dot{\theta}_2)\)が得られます。
\(
\begin{align}
L(r_1,\dot{r}_1,\theta_1,\dot{\theta}_1,& r_2,\dot{r}_2,\theta_2,\dot{\theta}_2) \\
=&\frac{1}{2}m_1(\dot{r}_1^2+r_1^2\dot{\theta}_1^2)+\frac{1}{2}m_2\left[\dot{r}_1^2+\dot{r}_2^2+r_1^2\dot{\theta}_1^2+r_2^2\dot{\theta}_2^2 \right. \\
&\left.+2(\dot{r}_1 r_2 \dot{\theta}_2-r_1\dot{r}_2 \dot{\theta}_1)\sin{(\theta_1-\theta_2)}+2(\dot{r}_1 \dot{r}_2 +r_1 r_2 \dot{\theta}_1 \dot{\theta}_2)+\cos{(\theta_1-\theta_2)}\right] \\
&+m_1gr_1\cos{\theta_1}+m_2g(r_1\cos{\theta_1}+r_2\cos{\theta_2})
\end{align}
\)

僕は先ほど式を簡単にするために座標変換をする、といいました。しかし、新しい座標系におけるラグランジアンはどう見ても元のデカルト座標系のラグランジアンに比べて複雑です。この理由は物理的な意味から来ています。
振り子をつないでいる棒が伸び縮みしないとすると系の自由度は角度\(\theta_1,\theta_2\)の２つです。
となると運動方程式は最高で2本の独立した方程式になるはずです。
しかし、デカルト座標の場合うまく自由度を落とすことができず、運動方程式は4つになってしまいます。
そこで棒が伸び縮みを簡単に表すことができる座標系に移ることでうまく方程式の数を減らせます。

新しい座標系でのラグランジアンで棒の伸び縮みがないという条件を表すには
\(
\begin{align}
r_1&=l_1\ (l_1\mbox{は定数}) \\
r_2&=l_2\ (l_2\mbox{は定数})
\end{align}
\)
と書けるわけで、また、
\(
\begin{align}
\dot{r}_1&=0 \\
\dot{r}_2&=0
\end{align}
\)
となるわけです。

\(m_1=m_2=m,\ l_1=l_2=l\)という場合を特に考えると、ラグランジアンは
\(
\displaystyle L(\theta_1,\dot{\theta}_1,\theta_2,\dot{\theta}_2)=ml^2\left[\dot{\theta}_1^2+\frac{1}{2}\dot{\theta}_2^2+\dot{\theta}_1\dot{\theta}_2\cos{(\theta_1-\theta_2)}\right]+mgl(2\cos{\theta_1}+\theta_2)
\)
と書けます。あとはラグランジュの運動方程式を当てはめて計算します。

3, 新たな座標系で運動方程式を導く
保存力場中でのラグランジュの運動方程式は
\(
\begin{eqnarray}
\left\{
\begin{aligned}
\frac{d}{dt}\left(\frac{\partial L}{\partial\dot{\theta}_1}\right)-\frac{\partial L}{\partial\theta_1}&=0 \\
\frac{d}{dt}\left(\frac{\partial L}{\partial\dot{\theta}_2}\right)-\frac{\partial L}{\partial\theta_2}&=0
\end{aligned}
\right.
\end{eqnarray}
\)
なので、代入し、\(\ddot{\theta}_1,\ddot{\theta}_2\)に関する運動方程式にすれば
\(
\begin{eqnarray}
\left\{
\begin{aligned}
\ddot{\theta}_1&=\frac{1}{2}\left\{-\ddot{\theta}_2\cos{(\theta_1-\theta_2)}-\dot{\theta}_2^2\sin{(\theta_1-\theta_2)}-2\frac{g}{l}\sin{\theta_1} \right\} \\
\ddot{\theta}_2&=-\ddot{\theta}_1\cos{(\theta_1-\theta_2)}+\dot{\theta}_1^2\sin{(\theta_1-\theta_2)}-\frac{g}{l}\sin{\theta_2}
\end{aligned}
\right.
\end{eqnarray}
\)
となります。これは非線形の２階連立微分方程式です。カオスです。解けません。
数値計算で解きます。用いるのは4次ルンゲ・クッタ法です。

実際に解いてみると、こうなります。

２重振り子の実験もされています。

[adsense2]

カオスが現れる場合、初期値鋭敏性という性質があります。
初期値がほんのちょっと変わるだけでその後の時間発展の様子が大きく変わる性質です。
ちょっと値を変えると上の動画とはまるで違う動きになります。
この動画は上の動画よりも中心に近いほうの粒子の初期速度が4%違います

機械精度の違いでさえも十分な時間発展の後には大きな違いが現れます。複雑であると同時に面白い現象です。

モンテカルロ法は，乱数を使って系の平衡状態にせまる一つの方法です。

ここでは二次元正方格子のイジング模型について紹介します。
イジング模型は隣り合う格子のスピン同士の相関エネルギーのみを考慮した模型です。
ハミルトニアンは以下のようになります。

\(
\displaystyle
H=-J \Sigma_{< i , j >}{s_{i}s_{j}}
\)
　
あるスピン状態を持つ系の状態Aについて，そのハミルトニアンを\(H_A\)とすると，
ボルツマン因子と分配関数\(Z\)を用いて状態の実現確率を計算することが出来ます。

\(
P(A)=\frac{\exp{(-\frac{H_A}{k_B T})}}{Z}
\)
　
Aからランダムに一つの格子におけるスピンを一つだけフリップした状態を状態Bとして，
状態Aと状態Bの実現確率比を乱数と比較します。
\(P(B)/P(A)>{\rm random}\)ならば状態Bを採用し，そうでなければ状態Aを採用します。
ここまでの流れを1mcs(モンテカルロステップ)といい，これを何度も繰り返すことで平衡状態に近づきます。
以下の二つのgifは，相転移温度より高い温度と低い温度で，
モンテカルロステップにより二次元正方格子のスピンがどのように振る舞うかをみているアニメーションです。
黒がダウンスピン，白がアップスピンです。

こちらは温度T=10eVで，自発磁化を持たない事が分かります。

こちらは温度T=0.001eVで，スピンが揃って自発磁化を持つ事が分かります。

画像をクリックするとmcsスタート!!

スピンのふるまいを追うと，自発磁化を持つ温度と持たない温度でスピンのふるまいが違う事が分かります。
二次元のイジング模型はオンサーガーにより厳密に解かれている(相転移温度T=2.26eV at J=1)ので，
それを踏まえて下のプログラムで遊んでみるといいですよ。

これのプログラムソースコード(保証は一切しません!!)
*注意　このプログラムを実行する前にising_anというディレクトリを作ってください。
*追記　2015 2/13　磁化とエネルギーを計算出来るようにしました。
計算方法は重み付き選択法というのを使ってます。

\(
\langle A \rangle = \frac{1}{T} \sum_{t=t_0+1}^{t_0 + T}A_{\alpha_{t}}
\)
　
十分平衡状態になっている\(t_0\)以降のモンテカルロステップで，物理量の統計平均を取る方法です。
平衡以前の物理量を計算しないmcsをからまわしと言います。

格子点51×51で，\(J=2\)で計算しました。
スピンの大きさを1として，相互作用をダブルカウントしているので注意しましょう。
即ち\(J=1\)の結果が欲しければプログラム上ではintJ=0.125にすればよいかも。

**********parameter**********
i_max::x方向の格子点の数
j_max::y方向の格子点の数
mcs_max::最大モンテカルロステップ数
mcs2step::mcs2step以降のmcsを平衡状態と見なし，物理量を計算します
intj::\(J/4\)
temp::温度
mov_max::mcs_max/movmax毎のスピン状態をgifアニメーションで見れるように区切ります。
**********parameter**********

単位格子あたりのエネルギーの温度依存性グラフを添付します。

相転移温度付近でエネルギーの変化が大きい事が分かります。

単位格子あたりの磁化の絶対値の温度依存性グラフを添付します。

相転移温度付近で自発磁化が発生し始めている事が分かります。

gnuplotを用いて，格子点上のスピンがどのように平衡状態へ近づくかをダイナミックに確認する事が出来ます。

gnuplot anime.gnuplot
　
を実行するとアニメが始まります。

gnuplot animegif.gnuplot
　
でgifアニメが作成出来ます。