この記事は、東北大学の「2022年度オープンキャンパス」のために作成された記事です。著作権は東北大学理学研究科物理学専攻 原子核理論研究室に帰属します。
原子核理論研究室のホームページに特設サイトを設置していますので、そちらも併せてご覧ください。
※数式を扱った説明が多いので、適宜飛ばしながらご覧ください。メインの数値計算は「3. 二重振り子とカオス」にあります。
目次
皆さんは「運動方程式」「ニュートン方程式」という言葉を聞いたことがありますか？高校で物理学を勉強している人なら一度は聞いたことがあるでしょう。高校物理を学んでいない人でも耳にしたことがあるかもしれません。
この「運動方程式」「ニュートン方程式」は、読んで字の如く、物体の運動を記述（説明）することのできる方程式で、1680年代にかの有名なイギリスの物理学者、アイザック・ニュートンによって作られた方程式です。この方程式のすごいところは、私たちの身の回りの物理（例えば投げた野球ボールの軌道）から天体の運動までを統一的に理解できるということです。
そんな運動方程式は次のような形をしています。\[F=ma\]この単純な1つの式が、さまざまな物理を説明することができるのです。
ではこの「運動方程式」を作っている\(m,a,F\)について具体的に見てみましょう。
まず、右辺に現れている\(m\)は、考えている物体の質量です。この量は秤（はかり）を用いれば簡単に知ることができますね。左辺にある\(F\)は、その物体に働いている力の大きさです。中学の理科の授業で力を表す単位として「ニュートン(N)」を習ったことがあると思います。その力を\(F\)と書きます。
では右辺にある\(a\)は何でしょうか。これは「加速度」と呼ばれる量で、「ある時間の間に、物体の速度がどれだけ変化しますか」という量です。例えば皆さんが歩いているときに、他の人に背中を押されたとしましょう。すると皆さんの歩く速度は少し早くなりますよね。つまり、力を加えられたことによって加速度が生じているのです。
もっと正確にいうと、加速度\(a\)は、ある時間\(\Delta t\)の間に速さが\(\Delta v\)だけ変化するとき
\[a=\frac{\Delta v}{\Delta t}\]
という式で与えられます。最初の速さが\(2\ \mathrm{m/s}\)で、\(0.1\)秒間背中を押されたことによって速さが\(3\,\mathrm{m/s}\)に変化したとしましょう。ここでの速さの変化は\(\Delta v=+1\ \mathrm{m/s}\)となります。時間の変化は\(\Delta t=0.1\ \mathrm{s}\)ですから、加速度は\(a=\Delta v/\Delta t = 10\ \mathrm{m/s^2}\)となります。
しかし、速度や加速度という量は目に見えるものではありません。座標（物体の位置）は目に見えてわかりやすいのに…そんな速度や加速度ですが、実は座標（と時間）を用いて書くことができます。これで加速度を使わずに運動方程式を書くことができますね！
まず、速度は「ある時間の間に、物体の座標がどれだけ変化するか」という量なので
\[v=\frac{\Delta x}{\Delta t}\]
とかけます。\(x\)というのは物体の座標で、一般に時間に依存して異なる値を持ちます。その意味で「\(x\)は時間\(t\)の関数」という見方もできるでしょう。
これらのことを用いると、運動方程式は次のように書くことができそうです。
\[m\frac{\Delta (\Delta x)}{(\Delta t)^2}=F\]
加速度\(a\)が少し気持ちの悪い形になってしまいました。これについて少し詳しくみてみましょう。\(\Delta x\)というのは\(x\)の変化量でした。つまり\(\Delta(\Delta x)\)は「\(x\)の変化量が、時間変化に伴ってどれだけ変化したか」という量に対応します。
これを具体的に考えてみます。ある時刻\(t_0\)において、物体が\(x_0\)という座標に存在しています。\(-\Delta t\)時間前の座標は\(x_-\)でした。この間の速度は
\[v_-=\frac{x_0-x_-}{\Delta t}\]
となります。また、\(\Delta t\)時間後の座標が\(x_+\)になるとしましょう。すると、ここでの速度は
\[v_+=\frac{x_+-x_0}{\Delta t}\]
です。この\(v_-,v_+\)という速度は、それぞれ時刻\(t_0-\Delta t/2\,,t_0+\Delta t/2\)の速度に対応します（\(t_0\pm \Delta t\)ではないことに気を付けてください）。このことから、速度の変化、加速度は次のように書くことができます。
\[a=\frac{v_+-v_-}{\Delta t}=\frac{x_+-2x_0+x_-}{(\Delta t)^2}\]
これは時刻\(t_0\)での加速度です。したがって、運動方程式は
\[m\frac{x_+-2x_0+x_-}{(\Delta t)^2}=F(t_0)\]
と表されます（\(F(t_0)\)は時刻\(t_0\)での力の大きさという意味です）。これによって、加速度を用いて表されていた運動方程式\(ma=F\)を、時刻と座標を用いた形に書き換えることができました！
先ほどまでは時間\(\Delta t\)の間に座標や速度がどれぐらい変わるかを考えました。ですが、このような方法で求めた速度・加速度は、もちろん\(\Delta t\)の選び方に依存してしまいます。同じ物体の運動を見ているのに、人によって速度や加速度の値が変わったら困りますよね。
なので\(\Delta t\)をとても小さく（無限に小さく）とった時の速度・加速度を考えます。
すでに述べましたが、座標\(x\)は時刻\(t\)の関数と見ることができます。そのため、時刻\(t\)から\(t+\Delta t\)までの座標の変化を\(\Delta t\)で割った量、すなわち速度\(\tilde v\)は
\[\tilde v=\frac{x(t+\Delta t)-x(t)}{\Delta t}\]
とかけます。今は\(\Delta t\)を\(0\)に近づけたいので、\(\Delta t\to 0\)という極限を取りましょう。すると、これは「微分」そのものです！
\[\lim_{\Delta t\to 0}\frac{x(t+\Delta t)-x(t)}{\Delta t}=\frac{dx(t)}{dt}=v(t)\]
これが「速度」（高校物理では「瞬間の速度」とも言います）の定義です。
同じように考えてみると、（瞬間の）加速度は速度の時間微分になります。
\[a(t)=\frac{dv(t)}{dt}=\frac{d^2x(t)}{dt^2}\]
このことから、運動方程式は
\[m\frac{d^2x(t)}{dt^2}=F(t)\]
という方程式になります。このように微分が含まれている方程式のことを「微分方程式」と言います。今の場合では、力の関数\(F(t)\)がわかっているときに、物体の軌道を表す関数\(x(t)\)がこの方程式を満たしている、関数\(x(t)\)は何ですか？という問題になるのです。
しかし一般的に、この微分方程式を正確に解く（\(x(t)\)を求める）ことはできません（できるケースはすこしだけあります）。つまり、力の関数がわかっても、運動方程式を眺めているだけでは物体がどのような運動をするのか知ることはできないのです！
そのために「数値計算」という技術が非常に重要になってくるのです。
では、運動方程式を解いて物体の運動を知るためにはどうすればよいでしょうか。そのために、すでに求めた次のような運動方程式の形を用います。
\[m\frac{x_+-2x_0+x_-}{(\Delta t)^2}=F(t_0)\]
ここで\(x_0\)はある時刻\(t_0\)における物体の座標で、\(x_+\,,x_-\)はそれぞれ時刻\(t_0+\Delta t\,,t_0 -\Delta t\)での座標でした。さてここで運動方程式を次のように変形してみましょう。
\[x_+=(\Delta t)^2\frac{F(t_0)}{m}+2x_0-x_-\]
この式が意味していることは何でしょうか？それは、「時刻\(t_0-\Delta t\)と\(t_0\)での座標、時刻\(t_0\)において物体にはたらいている力の大きさ\(F(t_0)\)がわかっていれば、右辺を計算することによって、時刻\(t_0+\Delta t\)における物体の座標\(x_+\)が決まる」ということです。ここで得られた\(x_+\)をさらに用いることで、さらに先の時刻\(t_0+2\Delta t\)での座標も求めることができます。
このようにして、「最初の座標」「その次の時刻での座標」の2つが分かれば、運動方程式から物体の座標を次々と求めることができるのです。
ここでは具体例として、物体の自由落下を考えてみましょう。地球上で自由落下する物体にはたらく力は
\[F=-mg\]
と書くことができます（符号が\(-\)なのは、力が下向きにはたらいているからです）。ここで\(g\)は重力定数\(g=9.8\)です。運動方程式は
\[x_+=-(\Delta t)^2g+2x_0-x_-\]
となりますので、これを用いて数値計算で物体の運動を計算します。
東京スカイツリー（高さ634m）から、初速度0で物体を落とす際は、
\[x_0=634, x_1=634\]
となります。\(x_0\)は最初の座標、\(x_1\)は\(\Delta t\)だけ時間が進んだ時の座標です。初速度を与えないという条件から、\(x_0=x_1=634\)としています。
次に載せているグラフは、時間幅\(\Delta t=0.1\ \mathrm{s}\)として、自由落下を計算した結果です。自由落下の運動の公式
\[x(t)=h-\frac{1}{2}gt^2\]
のグラフも載せていますので、数値計算で求められた結果と比較してみてください。
グラフ上の"Numerical"(赤色)は数値計算を用いた結果、"Exact"(青色)は自由落下の公式から得られる結果を示しています。赤色の線は青色の線とほとんど一致しており、数値計算を用いて運動の様子を比較的精度良く計算することができたことがわかります。
なお、今回はスカイツリーから物体を落とすケースを考えましたが、物体が地面に着いた後の跳ね返りなどは考慮していません。そのため、左のグラフでは時刻\(11\)秒後以降のデータが\(0\mathrm{[m]}\)になっており、右のグラフでは速度が\(-110\mathrm{[m/s]}\)から急激に\(0\mathrm{[m/s]}\)に変わっています。
1章では、私たちの身の回りで起こっている物理は、「運動方程式」という1つの式で説明できることを見ました。しかし、比較的単純な形をしている運動方程式で記述されるはずの私たちの身の回りの世界は、そう簡単に理解することができないほど複雑な現象で満ち溢れています。これはどういうことなのでしょうか。
その大きな理由の一つは、私たちの世界が「複数の物体が集まってできている」からです。私たちの身の回りに溢れているさまざまな「物質」は、すべて100数種類の「原子」からできています。そしてその「原子」は、「原子核」とその周りを回る「電子」からできています。「原子核」もまた、いくつかの「陽子」「中性子」からできていますし、「陽子」「中性子」は「クォーク」と呼ばれる「素粒子」からできています。つまり、私たちが普段目にしている「物質」は、膨大な量の小さな粒子が集まってできているのです。（「素粒子」などについての詳しい説明は、こちらの記事をご覧ください！）
一つ一つの粒子を記述する理論は単純でも、多くの粒子が集まると非常に複雑で多様な現象が現れます。これが、物理法則が単純な形をしているのに、私たちの身の回りにある現象が複雑である理由です。
1つの粒子を説明する理論を知っていても、粒子が複数集まった時に生まれる性質を直ちに理解することはできません。このような問題を「多体問題」と言います。私たち原子核理論研究室は、比較的単純な理論から生まれる、多体系の複雑な現象を理解することを目標に、日々研究をおこなっています。
ここでは多体問題をより身近に感じてもらうために、2つのトピックについてお話しします。
ここからは振り子が2つ繋がった「二重振り子」（下図）の運動方程式を計算して、初期状態が僅かに違うだけで最終的な結果が大きく変わってしまう現象（カオス）を見てみましょう。
Pendule_double.gif（著作者：Robert FERREOL パブリック・ドメイン）
二重振り子を考える前に、まずは単振り子の運動について考えてみましょう。振り子の振れ角を\(\theta\)とすると、運動方程式は
\[ml\frac{d^2\theta}{dt^2}=-mg\sin\theta\]
となります。高校物理では、\(\theta\)が十分に小さく\(\sin\theta\fallingdotseq\theta\)と近似できる場合の運動方程式
\[ml\frac{d^2\theta}{dt^2}=-mg\theta\]
を扱うことがほとんどでしょう。この場合には、運動方程式を解くと周期\(T=2\pi\sqrt{\frac{l}{g}}\)で振動する解（単振動）が得られます。この周期は振り子の重さ、振れ幅には依存しません。このことを「振り子の等時性」と言います。
しかし、これはあくまでも「\(\theta\)が十分に小さい」という仮定のもとで成り立つ性質です。では、一般の\(\theta\)の場合はどうなるでしょうか。この計算は大学生でも難しいのですが、結果は次のようになります。
\[T=4\sqrt{\frac{l}{g}}\int^{\pi/2}_0d\phi\frac{1}{\sqrt{1-\sin^2\frac{\theta_0}{2}\sin^2\phi}}\]
なお、\(\theta_0\)は振り子の最大角です。余談ですが、ここで現れている積分は「第一種完全楕円積分」と呼ばれる積分です。
では、高校で習う振り子の周期\(T=2\pi\sqrt{l/g}\)と、厳密に計算した周期はどれほど一致しているのでしょうか。それを示したのが下の図です。
\(\theta_0\)が大きくなるにつれて、"approximation" (\(\theta_0\)が小さいと近似した時の結果)よりも"exact"（近似を用いてない場合の結果）の方が大きくなることがわかります。ですが差があると言っても、\(\theta_0=40^\circ\)の時の差は\(0.0733\)程度。近似で求めた周期はおおよそ\(T=2\)ですので、大きく異なっているわけではなさそうです。
ここからは本題であった二重振り子を、「1. はじめに 〜物体の運動を数値計算で解く？〜」で説明した方法で計算していきます（高精度の計算が必要なので、全く同じ方法ではありません）。
運動を計算するためには運動方程式が必要です。単振り子の運動方程式は簡単でしたが、二重振り子の運動方程式はどうでしょうか。なお、二重振り子には2つの変数\(\theta_1\,,\theta_2\)がありますので、運動方程式は2つ出てきます。
その形はというと…
\[(m_{1}+m_{2})l_{1}\frac{d^2\theta_{1}}{dt^2}+m_{2}l_{2}\frac{d^2\theta_{2}}{dt^2}\cos(\theta _{1}-\theta _{2})+m_{2}l_{2}\left(\frac{d\theta_{2}}{dt}\right)^{2}\sin(\theta _{1}-\theta _{2})+(m_{1}+m_{2})g\sin \theta _{1}=0\]
\[ l_{1}l_{2}\frac{d^2\theta_{1}}{dt^2}\cos(\theta _{1}-\theta _{2})+l_{2}^{2}\frac{d^2\theta_{2}}{dt^2}-l_{1}l_{2}\left(\frac{d\theta_{1}}{dt}\right)^{2}\sin(\theta _{1}-\theta _{2})+gl_{2}\sin \theta _{2}=0\]
難しすぎる！この方程式を直接解くことは不可能です（数学的にも）。振り子の数が2つに増えただけで、運動を理解するのは難しくなってしまうのです。
ですが、数値計算を用いて運動を計算することはできます。次に載せているものは、複雑な運動方程式を数値的に解いて得られた振り子の運動です。初期条件である振り子の角度を色々と変えて計算してみてください。1つの振り子の初期の角度は、もう一つの振り子の角度よりも1度だけ大きく設定してあります。わずかな初期条件の違いが時間経過とともに大きくなる「カオス」の現象が見られるかもしれません。
※ここでは都合上、\(l_1=l_2=1,m_1=m_2=1\)として計算しています。