#0. はじめに
この記事はU-Tokyo mech (東京大学機械系) Advent Calendar 2018 25日目の記事です。
正確にはもう1月ですが、空いていたのを見つけたので、せっかくなので記事を書いておきました。
以前このツイートが少しバズっていたので、中身について解説しておこうと思います。
先日のカオス振り子動画を見て極端な場合を試してみたくなったので、3重振り子で高精度物理演算してみた。初期の振り子の角度を10の100乗(グーゴル)分の1ずつだけずらしたものを100個同時に再生。完全に一致してるように見えるけど、1分20秒くらいから「それ」がやってくる。 pic.twitter.com/a7E2OASIDe
— So Takamoto (JK) (@tkmtSo) 2017年5月25日
剛体重振り子はカオスが姿を現す例としてよく紹介されます。振り子という一見単純な動きに見えるものが予測不能な動きをする意外性、モデルの簡単さ、可視化の容易さ、何より現実で触って理解できるという点が魅力なのではないでしょうか。
この記事では、現象のカオス性については一旦脇に置いておいて、素朴な物理演算の例としてモデル化からコードを書くまでの流れを説明しようと思います。
読む人の対象としては、大学に入って力学を勉強した程度の人で、何かに応用をしてみたいなと思っている人を念頭に書いています。
#1. 方程式をつくる
まずは現象を記述する支配方程式を作ります。
今回目指すものは剛体3重振り子の運動です。剛体のリンク機構の古典的な運動をモデル化するので、ニュートンの運動方程式に基づいて支配方程式を作ればよいでしょう。
以降は振り子$i$の長さを$l_i$, 質量を$m_i$としておきます。
##きちんと解く
ところで、振り子なので隣り合った剛体はヒンジで点拘束されています。平面なので剛体の自由度は本来はx, y, 回転の3自由度ありますが、点拘束1つごとにxとy、つまり2つ分の自由度が落ちます。なので剛体3つ+ヒンジ3つの今回の系の自由度は実際には3つ(速度を入れて6つ)しかありません。
ニュートンの運動方程式を適用するには、自由度がどういう形であれ一旦デカルト座標上の関係に直して力のやりとりの関係から式を立てますが、この作業は結構面倒です。今回のように自由度が落ちた形で系を記述できる場合にはニュートンの運動方程式よりラグランジュの運動方程式を組み立てたほうがずっと楽になります。ラグランジュの運動方程式は位置エネルギーと運動エネルギーさえ書ければ自由度のとり方によらず運動方程式を立てることができます。
素直にパラメータをとるなら各振り子の角度 $\theta_1, \theta_2, \theta_3$ をパラメータとすればよいでしょう。各剛体棒の重心の位置$(x_i, y_i)$は、
\begin{align}
(x_1, y_1)&=\left(\frac{l_1}{2}\sin{\theta_1}, \frac{l_1}{2}\cos{\theta_1}\right) \\
(x_2, y_2)&=\left(l_1\sin{\theta_1}+\frac{l_2}{2}\sin{\theta_2}, l_1\cos{\theta_1}+\frac{l_2}{2}\cos{\theta_2}\right) \\
(x_3, y_3)&=\left(l_1\sin{\theta_1}+l_2\sin{\theta_2}+\frac{l_3}{2}\sin{\theta_3}, l_1\cos{\theta_1}+l_2\cos{\theta_2}+\frac{l_3}{2}\cos{\theta_3}\right)
\end{align}
剛体iの運動エネルギー$K_i$と位置エネルギー$U_i$は、質量$m_i$と慣性モーメント$I_i$を使って、
\begin{align}
K_i&=\frac{m_i}{2}\left(\dot{x_i}^2+\dot{y_i}^2\right)+\frac{I_i}{2}\dot{\theta_i} \\
U_i&=-m_igy_i
\end{align}
と書けます。$g$は重力加速度です。
また、棒が常識的な形(細くて一様)をしていれば慣性モーメントは $I_i=\frac{1}{12}m_il_i^2$ と書けます。
あとはラグランジュの運動方程式に従って$\theta_i$についての微分方程式を立てるだけです。機械系の人であれば大学の講義で2,3回遭遇しているとは思いますが、ラグランジュの運動方程式を以下に書いておきます。ある一般化座標$\boldsymbol{q}$とその時間微分$\dot{\boldsymbol{q}}$によって系の位置エネルギー$U$, 運動エネルギー$K$が書けるとき、ラグランジアン$L$は
L(\boldsymbol{q},\dot{\boldsymbol{q}})=K(\boldsymbol{q}, \dot{\boldsymbol{q}})-U(\boldsymbol{q}, \dot{\boldsymbol{q}})
と定義され、このときの運動方程式は、$\boldsymbol{q}$のi番目の要素を$q_i$として、
\frac{\partial}{\partial t}\left(\frac{\partial L}{\partial \dot{q_i}}\right)=\frac{\partial L}{\partial q_i}
と書けます。
さて、真面目に書き下してみましょう。
\begin{align}
L&=\left(\frac{1}{6}m_1+\frac{1}{2}m_2+\frac{1}{2}m_3\right)l_1^2\dot{\theta}^2+\left(\frac{1}{6}m_2+\frac{1}{2}m_3\right)l_2^2\dot{\theta_2}^2+\frac{1}{6}m_3l_3^2\dot{\theta_3}^2 \\
&+\left(\frac{1}{2}m_2+m_3\right)l_1l_2\cos{\left(\theta_1-\theta_2\right)}\dot{\theta_1}\dot{\theta_2}+\frac{1}{2}m_3l_1l_3\cos{\left(\theta_1-\theta_3\right)}\dot{\theta_1}\dot{\theta_3}+\frac{1}{2}m_3l_2l_3\cos{\left(\theta_2-\theta_3\right)}\dot{\theta_2}\dot{\theta_3} \\
&+g\left\{\left(\frac{1}{2}m_1+m_2+m_3\right)l_1\cos{\theta_1}+\left(\frac{1}{2}m_2+m_3\right)l_2\cos{\theta_2}+\frac{1}{2}m_3l_3\cos{\theta_3}\right\}
\end{align}
係数が煩雑なので、適当に記号を変換して書きなおすことにします。
\begin{align}
L=& \frac{1}{2}\left(\begin{array}{ccc}\dot{\theta_1} & \dot{\theta_2} & \dot{\theta_3}\end{array}\right)\left(\begin{array}{ccc}a_{11} & a_{12}\cos{\left(\theta_1-\theta_2\right)} & a_{13}\cos{\left(\theta_1-\theta_3\right)} \\ a_{12}\cos{\left(\theta_1-\theta_2\right)} & a_{22} & a_{23}\cos{\left(\theta_2-\theta_3\right)} \\ a_{13}\cos{\left(\theta_1-\theta_3\right)} & a_{23}\cos{\left(\theta_2-\theta_3\right)} &a_{33} \end{array}\right)\left(\begin{array}{c}\dot{\theta_1} \\ \dot{\theta_2} \\ \dot{\theta_3} \end{array}\right) \\
& + \left(\begin{array}{ccc}b_1 & b_2 & b_3 \end{array}\right)\left(\begin{array}{c}\cos{\theta_1} \\ \cos{\theta_2} \\ \cos{\theta_3} \end{array}\right)
\end{align}
このとき、ラグランジュの運動方程式は、
\left(\begin{array}{ccc}a_{11} & a_{12}\cos{\left(\theta_1-\theta_2\right)} & a_{13}\cos{\left(\theta_1-\theta_3\right)} \\ a_{12}\cos{\left(\theta_1-\theta_2\right)} & a_{22} & a_{23}\cos{\left(\theta_2-\theta_3\right)} \\ a_{13}\cos{\left(\theta_1-\theta_3\right)} & a_{23}\cos{\left(\theta_2-\theta_3\right)} & a_{33} \end{array}\right)\left(\begin{array}{c}\ddot{\theta_1} \\ \ddot{\theta_2} \\ \ddot{\theta_3} \end{array}\right) \\
=\left(\begin{array}{ccc}0 & -a_{12}\sin{\left(\theta_1-\theta_2\right)} & -a_{13}\sin{\left(\theta_1-\theta_3\right)} \\ a_{12}\sin{\left(\theta_1-\theta_2\right)} & 0 & -a_{23}\sin{\left(\theta_2-\theta_3\right)} \\ a_{13}\sin{\left(\theta_1-\theta_3\right)} & a_{23}\sin{\left(\theta_2-\theta_3\right)} & 0 \end{array}\right)\left(\begin{array}{c}\dot{\theta_1}^2 \\ \dot{\theta_2}^2 \\ \dot{\theta_3}^2 \end{array}\right)
+\left(\begin{array}{c}-b_1\sin{\theta_1} \\ -b_2\sin{\theta_2} \\ -b_3\sin{\theta_3} \end{array}\right)
となります。ふえー大変。数式処理ソフトは強い味方になってくれます。
運動方程式が行列の形で混ざって出てくるのはちょっと独特かもしれません。もちろん左から逆行列をかけることで左辺をシンプルにすることができますが、数式のままの計算ではちょっとやりたくないですね。
数値計算で解くときは、行列の中の値はすべて決まっていて逆行列を直接計算できるのでそんなに手間ではありません。$3\times 3$と決まっているならガウスの消去法を使うまでもありません。
ここまで一旦コードに落としてみましょう。以下はC++の例です。
#include <cmath>
#include <iostream>
struct mass_struct{
constexpr mass_struct(const double m0, const double m1, const double m2, const double l0,
const double l1, const double l2, const double g)
: a00((1.0/3.0*m0+m1+m2)*l0*l0), a11((1.0/3.0*m1+m2)*l1*l1), a22((1.0/3.0*m2)*l2*l2),
a01((0.5*m1+m2)*l0*l1), a02((0.5*m2)*l0*l2), a12((0.5*m2)*l1*l2),
b0(g*(0.5*m0+m1+m2)*l0), b1(g*(0.5*m1+m2)*l1), b2(g*(0.5*m2)*l2)
{}
double a00,a11,a22;
double a01,a02,a12;
double b0,b1,b2;
};
void ddt(const double (&theta)[3], const double (&dtheta)[3], const mass_struct& ms, double ans[3]){
const double t01 = theta[0]-theta[1];
const double t02 = theta[0]-theta[2];
const double t12 = theta[1]-theta[2];
const double cos01 = std::cos(t01), sin01 = std::sin(t01);
const double cos02 = std::cos(t02), sin02 = std::sin(t02);
const double cos12 = std::cos(t12), sin12 = std::sin(t12);
const double a01_cos01 = ms.a01*cos01, a01_sin01 = ms.a01*sin01;
const double a02_cos02 = ms.a02*cos02, a02_sin02 = ms.a02*sin02;
const double a12_cos12 = ms.a12*cos12, a12_sin12 = ms.a12*sin12;
const double dt0sq = dtheta[0]*dtheta[0];
const double dt1sq = dtheta[1]*dtheta[1];
const double dt2sq = dtheta[2]*dtheta[2];
const double g0 = -ms.b0*std::sin(theta[0])-a01_sin01*dt1sq-a02_sin02*dt2sq;
const double g1 = -ms.b1*std::sin(theta[1])+a01_sin01*dt0sq-a12_sin12*dt2sq;
const double g2 = -ms.b2*std::sin(theta[2])+a02_sin02*dt0sq+a12_sin12*dt1sq;
//calc 3x3 inverse matrix
const double det_a_inv = 1.0/(ms.a00*ms.a11*ms.a22+2.0*a01_cos01*a02_cos02*a12_cos12
-ms.a00*a12_cos12*a12_cos12-ms.a11*a02_cos02*a02_cos02
-ms.a22*a01_cos01*a01_cos01);
const double a_inv_0 = ms.a11*ms.a22-a12_cos12*a12_cos12;
const double a_inv_1 = ms.a00*ms.a22-a02_cos02*a02_cos02;
const double a_inv_2 = ms.a00*ms.a11-a01_cos01*a01_cos01;
const double a_inv_01 = a02_cos02*a12_cos12-a01_cos01*ms.a22;
const double a_inv_02 = a01_cos01*a12_cos12-a02_cos02*ms.a11;
const double a_inv_12 = a02_cos02*a01_cos01-a12_cos12*ms.a00;
ans[0] = det_a_inv*(a_inv_0 *g0+a_inv_01*g1+a_inv_02*g2);
ans[1] = det_a_inv*(a_inv_01*g0+a_inv_1 *g1+a_inv_12*g2);
ans[2] = det_a_inv*(a_inv_02*g0+a_inv_12*g1+a_inv_2 *g2);
}
int main(int argc, char** argv){
double theta[3] = {1.0, 0.1, 0.1};
double dtheta[3] = {0.0, 0.0, 0.0};
constexpr double l0 = 1.0, l1 = 1.0, l2 = 0.5;
constexpr double m0 = 1.0, m1 = 1.0, m2 = 0.5;
constexpr double g = 1.0;
constexpr auto ms = mass_struct(m0, m1, m2, l0, l1, l2, g);
double ddtheta[3] = {};
ddt(theta, dtheta, ms, ddtheta);
std::cout << "ddt: " << ddtheta[0] << " " << ddtheta[1] << " " << ddtheta[2] << std::endl;
return 0;
}
長さが$1:1:0.5$になっているのは単に大学1年のときの物理実験室に置いてあった振り子のおもちゃの長さがこうだったからです。
添字が数式とずれてしまっているので注意してください。今回のようにざくざく計算をしようみたいな発想のときはC++で書きたくなりますが、べたっと書いてるのでちょっと見づらいかもしれません。
ともあれ微分方程式を得ることができました。
##おまけ: もう少しゆるっとモデル化する
なんでこんなに大変になったかというと、隣り合う振り子の拘束条件が幾何的にややこしい設定だったからでした。
そこで、完全に点拘束にする代わりに、接合部を「ものすごく硬いバネで結ぶ」モデルで置き換えることを考えてみます。
何が起きるかを箇条書きにすると、
- 点拘束は完全に守られるものではなく、運動の結果として「ほぼ」守られるものになる
- 系の自由度は3ではなく9になる
- ラグランジュの運動方程式を使わなくても各剛体にかかる力から素直にニュートンの運動方程式が立式できる
- ばねの硬さという新しい値を設定する必要がある
- エネルギー保存則は引き続き守られる
などがあります。
式で書くと以下のようになります。剛体棒の位置と回転を $(x_i, y_i, \theta_i)$ で表すと、i番目のばねの力$\boldsymbol{f_i}$は、ばね定数$k_i$を用いて
\begin{align}
f_{ix}&=k_i\left\{\left(x_i-\frac{1}{2}l_i\sin{\theta_i}\right)-\left(x_{i-1}+\frac{1}{2}l_{i-1}\sin{\theta_{i-1}}\right)\right\} \\
f_{iy}&=k_i\left\{\left(y_i-\frac{1}{2}l_i\cos{\theta_i}\right)-\left(y_{i-1}+\frac{1}{2}l_{i-1}\cos{\theta_{i-1}}\right)\right\}
\end{align}
ただし0番目の剛体棒は原点の固定点としています。
このばねがi番目の剛体棒にかける力は$\boldsymbol{f_i}$の$\theta_i$方向成分、モーメントは$\boldsymbol{f_i}$の$\theta_i$に垂直方向成分ということですぐ求まります。他の力は重心にかかる重力だけです。このまま各剛体ごとに運動方程式が立ちます。運動方程式は、剛体にかかる力とモーメントを$\boldsymbol{F_i}, M$として、
\begin{align}
m_i\ddot{x_i}&=F_{ix} \\
m_i\ddot{y_i}&=F_{iy} \\
I_i\ddot{\theta_i}&=M_i
\end{align}
シンプルですね。力のやりとりが局所的になったため、各要素ごとに独立に運動を解くことができるようになりました。これは逆行列が出てこなくなったこととも関係しています。運動方程式を導出する点ではこっちのモデル化のほうが楽でしょう。
ラグランジュ方程式を解かなくてよくなったというのは拡張性の点からもうれしいことがあります。振り子の数を変えても同じやり方で普通に解けるのです。
先のC++のコードがvectorなどを使わず剛体3個で決め撃ちして書いている理由の1つにもなっています。先のモデル化では運動方程式を剛体の数の増減に対応させるのが極めて難しいのです。
また、曖昧な表現ですが、数値解析としてはなんとなくロバストにもなります。動ける範囲が位相空間で3次元の薄い多様体だったのが9次元の厚みのある(?)多様体になった、といえばいいでしょうか。これは実用上も使われていて、例えば有限要素法における物体間の接触では、表面間の幾何的な接触判定条件を厳密に解く代わりに表面同士をバネで反発させる方法が使われることがあります。
計算上でも特性の違いが出てきます。次の章の最後に少し書いておきます。
#2. 離散化する
オイラー法とかルンゲクッタ法とかの話をします。(時間依存の話題は抜いています。)
##汎用的な処方箋
なんだかんだで微分方程式が立式できましたが、コンピュータはたかだか有限回の計算しかできないため、このままでは解けません。そこで微分方程式を差分方程式で近似する離散化という操作を行います。
とりあえず2階の微分方程式を1階の微分方程式にしましょう。今のところ、以下のような連立微分方程式が得られています。
\frac{\mathrm{d}^2}{{\mathrm{d}t}^2}x=f\left(\frac{\mathrm{d}x}{\mathrm{d}t}, x\right)
変数の数を増やしてよいならこれは簡単で、以下のようにします。
\frac{\mathrm{d}}{\mathrm{d}t}\left(\begin{array}{c}
\dot{x} \\
x
\end{array}\right) =
\left(\begin{array}{c}
f\left(\dot{x}, x\right) \\
\dot{x}
\end{array}\right)
$\dot{x}$と$x$の組を$x$と置き直して、結局、$\mathrm{d}x/\mathrm{d}t=f(x)$の形になります。
ここから先が問題です。$dt$を有限の幅の$\Delta t$に置き換えた式が欲しいのです。
おそらく最もシンプルな考え方は、上の式の$\mathrm{d}t$を直接$\Delta t$に置き換える方法です。おそらく以下のようになります。
\begin{align}
k_1&=f(x(t)) \\
x(t+\Delta t)&=x(t)+k_1\Delta t
\end{align}
これを有限小の間隔$\Delta t$ごとに繰り返します。このやり方は(陽的)オイラー法と呼ばれています。実装してみるとわかりますが、$\Delta t$をいくら小さくとっても精度が悪いです。
関数の微分値を求めるときに微小な$\Delta t$をとって数値的に求めるといい値が出るので、逆向きの試みがうまくいかないのは少し意外でもあります。
1階微分でうまくいかない場合はどうするかというと、当然2階微分を使うという発想になります。
2階微分を求めるには「変化の変化ぶん」が必要になります。つまり2点は計算が必要なことがわかります。傾きの計算を2回繰り返して平均をとる方法が思いつきます。
\begin{align}
k_1&=f(x(t)) \\
k_2&=f(x(t)+k_1\Delta t) \\
x(t+\Delta t)&=x(t)+\frac{1}{2}\left(k_1+k_2\right)\Delta t
\end{align}
これは2次のルンゲクッタ法あるいはホイン法といい、実際うまくいきます。直感的にはこれで2次関数的なカーブを追えるからです。
このやり方でもっと高階の微分まで使おうというのが、数値解析の分野で広く使われている(4次の)ルンゲクッタ法です。4点とって4次の項まで合わせにいきます。
複数の補間点をどの幅でとってどう混ぜればいいのかというのは結構難しい問題です。というのも係数の選び方に任意性があるからです。一般的には以下の(古典的な)ルンゲクッタ法の式が使われます。
\begin{align}
k_1&=f(x(t)) \\
k_2&=f(x(t)+\frac{1}{2}k_1\Delta t) \\
k_3&=f(x(t)+\frac{1}{2}k_2\Delta t) \\
k_4&=f(x(t)+k_3\Delta t) \\
x(t+\Delta t)&=x(t)+\frac{1}{6}\left(k_1+2k_2+2k_3+k_4\right)\Delta t
\end{align}
これでテイラー展開の4次の項まで一致させたカーブで微分方程式を解くことができます。
実装は以下のようになります。性能が良い悪いというのは数字で見たほうがはっきりします。試しに$\mathrm{d}x/\mathrm{d}t=x, x(0)=1$の微分方程式をオイラー法とルンゲクッタ法で比較してみましょう。厳密解は$x=\exp{(t)}$です。
#include <iostream>
#include <iomanip>
#include <cmath>
#include <utility>
double f(const double x){
return x;
}
double euler(const double x, const double dt){
return x + dt * f(x);
}
double runge_kutta(const double x, const double dt){
const double k1 = f(x);
const double k2 = f(x+0.5*dt*k1);
const double k3 = f(x+0.5*dt*k2);
const double k4 = f(x+dt*k3);
return x + dt / 6.0 * (k1 + 2.0*k2 + 2.0*k3 + k4);
}
auto compare(const int n_iter, const double final_t){
const double dt = final_t/n_iter;
double xe = 1.0;
double xr = 1.0;
for(int i = 0; i < n_iter; i++){
xe = euler(xe, dt);
xr = runge_kutta(xr, dt);
}
const double exact = std::exp(n_iter*dt);
return std::make_tuple(dt, xe, xr);
}
int main(){
const double final_t = 10.0;
const double exact = std::exp(final_t);
std::cout << "dx/dt=x, x(0)=1.0 -> x(" << final_t << ")=" << exact << std::endl;
for(int i = 0; i < 5; i++){
const int n_iter = 10 * std::pow(10, i);
double dt, euler, runge;
std::tie(dt, euler, runge) = compare(n_iter, final_t);
std::cout << std::setprecision(6) << "dt= " << std::setw(8)
<< dt << ", N= " << std::setw(8) << n_iter
<< ", Euler error: " << std::setw(13) << euler/exact - 1.0
<< ", Runge-Kutta error: " << std::setw(13) << runge/exact - 1.0 << std::endl;
}
return 0;
}
実行結果はこちら。
dx/dt=x, x(0)=1.0 -> x(10)=22026.5
dt= 1, N= 10, Euler error: -0.95351, Runge-Kutta error: -0.0360016
dt= 0.1, N= 100, Euler error: -0.374361, Runge-Kutta error: -7.66777e-06
dt= 0.01, N= 1000, Euler error: -0.0484558, Runge-Kutta error: -8.26415e-10
dt= 0.001, N= 10000, Euler error: -0.00498421, Runge-Kutta error: -8.09353e-14
dt= 0.0001, N= 100000, Euler error: -0.000499842, Runge-Kutta error: -1.38778e-14
オイラー法も変な値に漸近するわけではなく、10万分割すれば3ケタくらい合います。一方ルンゲクッタ法はたったの100分割でもその精度を上回っていて、1万分割もすると倍精度浮動小数点数の打切り誤差のオーダーに到達してしまいます。
微分値の計算の数はオイラー法のたったの4倍にしかなっていないことを考えると素晴らしい性能です。
ところで、5次以降のルンゲクッタ法ですが、おもしろいことに5次のルンゲクッタ法を使うには補間点が5点だけでは不可能なことが知られています。
##おまけ: 運動方程式向きの処方箋
実際の機械系の数値解析の分野で運動方程式を解くときにいつもルンゲクッタ法が使われるかというと、そうでもありません。
例えば連続体を扱う有限要素法ではNewmarkのβ法という陰的な(右辺に$x(t+\Delta t)$が出てくる)解法が使われることがあります。陰的な解法は前向きに解けないため実装の点からも計算コストの点からも大変ですが、収束性が良く大きな$\Delta t$をとっても破綻しないなどのメリットがあります。これは有限要素法ではうれしい性質です。
また分子動力学法では速度ベルレ法という、一見するとホイン法の亜種のような2次精度の手法をよく使います。計算手順としては先に速度を更新してから位置を更新するという程度のシンプルな方法ですが、速度ベルレ法はなんと系のエネルギーを長期的に保存するという顕著な特性があります。分子動力学法は統計集団を扱う都合上、系の軌跡の正確さそのものよりも系全体のマクロな量(エネルギーとか)が保存されることが重要なため、この性質は都合が良いのだと思います。
有限な$\Delta t$をとっているのに誤差が蓄積しないというのは不思議な性質です。速度ベルレ法はシンプレクティック法と呼ばれる手法の一種で、この特性にはハミルトン力学とシンプレクティック性に関する深遠な話題が裏にあります。興味がある人は調べてみてください。
ちなみに、オイラー法やルンゲクッタ法のような陽的解法で保存系をシミュレーションするとエネルギーが次第に上昇していきます。ではなぜ今回ルンゲクッタ法を使ったかというと、実はシンプレクティック法は適用可能な対象が限られている(ひらたく書くと、位置エネルギーが$q$だけ、運動エネルギーが$\dot{q}$だけで書かれている必要がある)からです。その他にも$\Delta t$を一定にしないといけなかったりと、独特の適用条件があります。科学技術の世界における「制約が強い高位の魔法」みたいなものの一つだと思っています。
(今回のような場合でも為す術がないわけではなく、正準変換してハミルトンの運動方程式にすれば適用可能になります。手間はかかりますが。)
実は「もう少しゆるっとモデル化する」のところで説明したモデルはこの条件を満たしています。つまり速度ベルレ法が適用可能だったりします。
こちらのツイートの動画はこのゆるっとモデル+速度ベルレ法の組み合わせに変更して運動方程式を解いています。こちらのほうがエネルギーが安定しています。
振り子の形と計算手法を手直しして再計算。カオス的な系では先の例のように未来が"爆発"するけど、条件をうまく選べばこのように無数の現在を物理演算を行った先の未来で収束させることができる。 pic.twitter.com/YxPKs1j2eo
— So Takamoto (JK) (@tkmtSo) 2017年5月27日
#3. 可視化する
可視化は物理シミュレーションの醍醐味です。しかし慣れないとハードルが高く感じるパートでもあると思います。今回のように特定の対象向けの可視化専用ツールにはまらないようなときは特にそう見えるでしょう。
ゲームのように実時間でシミュレーションする面白さもあるのですが、まずは解析結果をファイルに出力しておいて後でまとめて可視化したほうが楽です。
なんらかの空間分布なり配置を持つ物理シミュレーションの結果を可視化する汎用的なツールとして、ParaViewというソフトウェアがあるので紹介しておきます。
ParaViewでは、3次元空間に配置されるデータを読み込み、値に応じて3D形状を作成したり切り取ったり彩ったりといった加工をしていくことができます。
流体のようなメッシュ状のデータでも、頂点や直線やポリゴンからなるデータでも対応しています。データは様々な形式に対応していますが、例えば.vtkファイルを通せば素朴に読み書きができると思います。形式に従って順番に頂点、直線、ポリゴンのデータを書いていくだけです。
ParaView側で適当に色を付けたり頂点や直線に球体や柱を割り当てたりすれば見栄えがよくなります。
横から見たらどうなってるのっと pic.twitter.com/Rg1AwrSUN2
— So Takamoto (JK) (@tkmtSo) 2017年5月27日
フォトリアリスティックな描画がしたかったらParaViewで作成した3次元メッシュをさらに他のCGソフト(例えばBlenderなど)に流し込みましょう。
#4. 実装する
##ソースコード
ソースコードの全体はこんな感じになります。おまけ要素として$\theta$が$0$から$2\pi$の範囲に収まるように切り取っています。特になくても動きます。
このままコンパイルして直下にoutディレクトリを作って実行して、outディレクトリ以下のvtkファイルをParaViewで読ませたらとりあえず動画になります。
(そこそこ古くないC++のコンパイラが必要になります。例えばGCCならC++11に対応するように-std=c++11オプションをつけてコンパイルしましょう。)
もちろん、3重振り子に限らずもっと汎用的なソルバにしたければ要素数を3に限定しないように書き直したほうがいいでしょう。
#include <cmath>
#include <string>
#include <iostream>
#include <fstream>
const double PI = std::atan(1.0)*4;
struct mass_struct{
constexpr mass_struct(const double m0, const double m1, const double m2,
const double l0, const double l1, const double l2, const double g)
: a00((1.0/3.0*m0+m1+m2)*l0*l0), a11((1.0/3.0*m1+m2)*l1*l1), a22((1.0/3.0*m2)*l2*l2),
a01((0.5*m1+m2)*l0*l1), a02((0.5*m2)*l0*l2), a12((0.5*m2)*l1*l2),
b0(g*(0.5*m0+m1+m2)*l0), b1(g*(0.5*m1+m2)*l1), b2(g*(0.5*m2)*l2)
{}
double a00,a11,a22;
double a01,a02,a12;
double b0,b1,b2;
};
void ddt(const double (&theta)[3], const double (&dtheta)[3], const mass_struct& ms, double ans[3]){
const double t01 = theta[0]-theta[1];
const double t02 = theta[0]-theta[2];
const double t12 = theta[1]-theta[2];
const double cos01 = std::cos(t01), sin01 = std::sin(t01);
const double cos02 = std::cos(t02), sin02 = std::sin(t02);
const double cos12 = std::cos(t12), sin12 = std::sin(t12);
const double a01_cos01 = ms.a01*cos01, a01_sin01 = ms.a01*sin01;
const double a02_cos02 = ms.a02*cos02, a02_sin02 = ms.a02*sin02;
const double a12_cos12 = ms.a12*cos12, a12_sin12 = ms.a12*sin12;
const double dt0sq = dtheta[0]*dtheta[0];
const double dt1sq = dtheta[1]*dtheta[1];
const double dt2sq = dtheta[2]*dtheta[2];
const double g0 = -ms.b0*std::sin(theta[0])-a01_sin01*dt1sq-a02_sin02*dt2sq;
const double g1 = -ms.b1*std::sin(theta[1])+a01_sin01*dt0sq-a12_sin12*dt2sq;
const double g2 = -ms.b2*std::sin(theta[2])+a02_sin02*dt0sq+a12_sin12*dt1sq;
//calc 3x3 inverse matrix
const double det_a_inv = 1.0/(ms.a00*ms.a11*ms.a22+2.0*a01_cos01*a02_cos02*a12_cos12
-ms.a00*a12_cos12*a12_cos12-ms.a11*a02_cos02*a02_cos02
-ms.a22*a01_cos01*a01_cos01);
const double a_inv_0 = ms.a11*ms.a22-a12_cos12*a12_cos12;
const double a_inv_1 = ms.a00*ms.a22-a02_cos02*a02_cos02;
const double a_inv_2 = ms.a00*ms.a11-a01_cos01*a01_cos01;
const double a_inv_01 = a02_cos02*a12_cos12-a01_cos01*ms.a22;
const double a_inv_02 = a01_cos01*a12_cos12-a02_cos02*ms.a11;
const double a_inv_12 = a02_cos02*a01_cos01-a12_cos12*ms.a00;
ans[0] = det_a_inv*(a_inv_0 *g0+a_inv_01*g1+a_inv_02*g2);
ans[1] = det_a_inv*(a_inv_01*g0+a_inv_1 *g1+a_inv_12*g2);
ans[2] = det_a_inv*(a_inv_02*g0+a_inv_12*g1+a_inv_2 *g2);
}
void runge_kutta(double (&theta)[3], double (&dtheta)[3], const mass_struct& ms, const double dt){
double dk1[3];
ddt(theta, dtheta, ms, dk1);
// double k2[3] = {dtheta[0], dtheta[1], dtheta[2]}; // same as dtheta
double tmp[3], dtmp[3];
tmp[0] = theta[0]+0.5*dt*dtheta[0];
tmp[1] = theta[1]+0.5*dt*dtheta[1];
tmp[2] = theta[2]+0.5*dt*dtheta[2];
dtmp[0] = dtheta[0]+0.5*dt*dk1[0];
dtmp[1] = dtheta[1]+0.5*dt*dk1[1];
dtmp[2] = dtheta[2]+0.5*dt*dk1[2];
double dk2[3];
ddt(tmp, dtmp, ms, dk2);
// double k2[3] = {dtheta[0], dtheta[1], dtheta[2]}; // same as dtheta
// tmp[0] = theta[0]+0.5*dt*k2[0]; // do not change
// tmp[1] = theta[1]+0.5*dt*k2[1];
// tmp[2] = theta[2]+0.5*dt*k2[2];
dtmp[0] = dtheta[0]+0.5*dt*dk2[0];
dtmp[1] = dtheta[1]+0.5*dt*dk2[1];
dtmp[2] = dtheta[2]+0.5*dt*dk2[2];
double dk3[3];
ddt(tmp, dtmp, ms, dk3);
// double k3[3] = {dtheta[0], dtheta[1], dtheta[2]}; // same as dtheta
tmp[0] = theta[0]+dt*dtheta[0];
tmp[1] = theta[1]+dt*dtheta[1];
tmp[2] = theta[2]+dt*dtheta[2];
dtmp[0] = dtheta[0]+dt*dk3[0];
dtmp[1] = dtheta[1]+dt*dk3[1];
dtmp[2] = dtheta[2]+dt*dk3[2];
double dk4[3];
ddt(tmp, dtmp, ms, dk4);
// double k4[3] = {dtheta[0], dtheta[1], dtheta[2]}; // same as dtheta
theta[0] += dt*dtheta[0];
theta[1] += dt*dtheta[1];
theta[2] += dt*dtheta[2];
dtheta[0] += dt/6.0*(dk1[0] + 2.0*dk2[0] + 2.0*dk3[0] + dk4[0]);
dtheta[1] += dt/6.0*(dk1[1] + 2.0*dk2[1] + 2.0*dk3[1] + dk4[1]);
dtheta[2] += dt/6.0*(dk1[2] + 2.0*dk2[2] + 2.0*dk3[2] + dk4[2]);
if(theta[0] < 0.0) theta[0] += 2.0*PI;
else if(theta[0] > 2.0*PI) theta[0] -= 2.0*PI;
if(theta[1] < 0.0) theta[1] += 2.0*PI;
else if(theta[1] > 2.0*PI) theta[2] -= 2.0*PI;
if(theta[2] < 0.0) theta[2] += 2.0*PI;
else if(theta[2] > 2.0*PI) theta[2] -= 2.0*PI;
}
int main(int argc, char** argv){
double theta[3] = {1.6, 1.6, 1.6};
double dtheta[3] = {0.0, 0.0, 0.0};
constexpr double l0 = 1.0, l1 = 1.0, l2 = 0.5;
constexpr double m0 = 1.0, m1 = 1.0, m2 = 0.5;
constexpr double g = 1.0;
constexpr auto ms = mass_struct(m0, m1, m2, l0, l1, l2, g);
const int out_interval = 100;
for(int i = 0; i < 200000; i++){
runge_kutta(theta, dtheta, ms, 0.001);
if(i%out_interval == 0){
std::string outstr = "out/test."+std::to_string(i/out_interval)+".vtk";
std::fstream fw(outstr.c_str(), std::ios::out);
fw << "# vtk DataFile Version 2.0\nTest1\nASCII\nDATASET POLYDATA\nPOINTS 4 float\n";
fw << 0.0 << " " << 0.0 << " 0.0\n";
fw << l0*std::sin(theta[0]) << " " << l0*std::cos(theta[0]) << " 0.0\n";
fw << l0*std::sin(theta[0])+l1*std::sin(theta[1]) << " "
<< l0*std::cos(theta[0])+l1*std::cos(theta[1]) << " 0.0\n";
fw << l0*std::sin(theta[0])+l1*std::sin(theta[1])+l2*std::sin(theta[2]) << " "
<< l0*std::cos(theta[0])+l1*std::cos(theta[1])+l2*std::cos(theta[2]) << " 0.0\n";
fw << "LINES 3 9\n";
fw << "2 0 1\n";
fw << "2 1 2\n";
fw << "2 2 3";
fw.close();
}
}
return 0;
}
##精度の話
今回のコードはdouble型(64bit浮動小数点型)で計算しています。はじめのツイートでやっていたのは、この中身の計算をすべて任意精度演算に置き換えて計算したものです。(GMP: GNU Multi-precision Libraryを使って512bitの精度で計算しています。)
今回のように四則演算をほぐして書いておくと、可読性という点では微妙ですが、GMPへの置き換えという点では便利かと思います。
数値計算上の誤差という点では誤差を極めて小さくできているため、カオスに特有の初期値鋭敏性(はじめのごく微小な差が時間とともに拡大し続ける)が再現できています。特に、$10^{-100}$という極めて小さいオーダーの初期値のずれを入れられたことで、肉眼レベルの分岐をずっと後になってから起こすことに成功しています。
ではこの計算は合っているのかというと、そんなことはありません。空間方向(値の大小)こそ精度を上げていますが、離散化の段階で時間方向の$\Delta t$が有限の小さくない値になっています。このため、得られた解はとても微分方程式の真の解と言えるものにはなっていません。(先のツイートでは高精度とうたってしまったので、書いた後でしまったこれはミスリーディングな表現になってしまったなあと思っていました。)
このあたりをきちんと扱いたければ、計算精度をきちんと保証しながら数値積分するなどのより高度な手法に頼る必要があるでしょう。逆に、振り子のモデルとしてはむしろ予測が実質的に不可能であるということこそが面白い点かなと思います。
他に精度にまつわる話としては、現実の振り子は剛体ではないという話があります。ここまで来るとなんだか揚げ足を取るような話ですが、例えばこのモデルでは振り子の端にかかった力が(光の速さを超えて)一瞬で振り子の反対側まで届くため、何かしらを近似したモデルなのだなということはわかるかと思います。
#5. まとめ
こういったさくっとシミュレーションが手軽にできるようになると研究も遊びも幅が広がっていいと思います。興味があったらぜひ遊んでみてください。
24日のアドベントカレンダーは@kn1chtさんのLaTeXで書いた学位論文もtextlintでチェックしよう（&工学系論文向けのtextlint presetを作っている話）です。