この記事は KLab Engineer Advent Calendar 2021 24日目の記事です。
Qiitaへの投稿は一昨年の12/24、去年の12/24のアドカレに続いて3回目になります suzuna-honda といいます。よろしくお願いします。
はじめに
昨今はゲームエンジン側に標準で物理シミュレーションが搭載されており、ゲーム開発中に自前で物理挙動を実装することは少なくなってきたかと思います。とはいえ、では物理学について何もしらなくて良いかというと、決してそうではないと私は考えます。最低限の基礎知識はあってこそ、物理シミュレーションの勘所を知り有効活用ができるのではないでしょうか。
というわけで、簡単なところで古典力学から改めて勉強していきましょう。
（以下、元々Qiitaエントリ用に作成した文章ではないため、図が存在しない/怪しい点など多々あるかと思いますが、優しい心でご覧ください）
chapter 1 質点の運動力学.
本章では、古典力学の中でも基礎中の基礎となる"質点"を用いて物体の運動を考えていきます。質点は後に紹介する"剛体"の素地でもあり、質点をベースとした物理シミュレーションもありますので、最後までお世話になる重要な概念です。この質点を用いて、運動とは何なのか、どう表現するのか、そして何故"ニュートンの運動法則"のような基本法則が生まれるに至ったのか、その逐一を順次解説していきます。
1.1 質点の運動学.
1.1.1 質点とは.
物事をシンプルに捉えるために、まずは"質点"というモデルを導入します。質点とは、何らかの物体を、大きさを無視した"点"であると見立てて、その点に物体の全部の重さが集まっている、とした仮想的な物体です。質点は向き/姿勢というパラメータを持たないので、並進運動のみが行われます。
例えば野球の試合中、バッターの打った打球がバットから離れた瞬間からどのような軌跡を描くか、また打球が野球場の何処に落ちるかを知りたいとします。このとき、打球自体の大きさは無視して点と見立てて、その点がどのように動くのか扱うようにしても差し支えはなさそうです。このときの打球のように、大きさを考慮に入れない物体の運動を考える場合、"質点"を用いるわけです。
［ 図1.1 ］野球ボールの中心に点を描いて、その点＝質点の軌跡を放物運動で記した絵。
むろん、実際の打球には大きさがあり、向き/姿勢があり、回転しています。この回転によって、よく言われるような「打球が伸びる」などといったことが起こったりします。が、そういった煩雑な事象に最初から立ち向かっては身が持ちません。複雑な現実世界を理解するためにも、まずは簡単なことから積み重ねていくべきです。このように物事を単純化して考えることを"理想化"といいます。
1.1.2 位置と時間.
質点が今どこにあるのかを一意に表すため、まず適当な座標系を用意し、その上で位置ベクトル＜ｒ = （ x,y,z ）＞を用いて表現します。
［ 図1.2 ］座標軸x3と質点の位置、そこに伸びる位置ベクトルの絵。
質点が動いているならば、時間が経つにつれてこの＜ｒ＞の値は変化していきます。その度々に"時刻＜t＞での位置＜ｒ＞"と記述していては分かりづらく不便です。そこで、指定した時刻＜t＞での質点の位置＜ｒ＞を＜ｒ＝ｆ（t）＞、短縮して＜ｒ（t）＞という関数の形で記述してやることで、一意な値として定めることが出来ます。つまり質点の"いつ"の位置ベクトルが欲しいのか、時刻という変数＜t＞を渡してやってはじめて＜ｒ＞の値が確定するわけです。＜ｒ（t）＞はベクトルですので、それぞれの各成分は＜ｘ（t）,ｙ（t）,ｚ（t）＞となります。この時間の関数＜ｒ（t）＞を全ての時刻毎に求めていくと、その質点の運動の軌跡が分かります。
［ 図1.3 ］図1.2の質点を動かして、いくつかの点が並んでそれが曲線となり、ｒ（0-n）となる絵。
＜ｒ（t）＞と毎回律儀に記述していると肝心の式が見づらくなってしまう場合がありますので、時刻＜t＞が流れから読みとれる場合には式から省略しても構わないとします。その場合であっても、位置ベクトル＜ｒ＞が時間の関数であるという事実は変わりません。
1.1.3 速度.
図1.4のように、坂の下に砂地のある滑り台からボールを転がす様子を観察してみましょう。ボールを頂上に置き、軽く蹴っ飛ばしてやります。最初はそのままゴロゴロと転がっていきますが、下り坂で加速をはじめ、坂が終わったところで砂地に到達し、砂に埋まったのか急激に減速してボールは止まりました。初期位置から止まった場所までの距離が10メートル、止まるまでに5秒掛かりました。このボールの、転がっていた時の速さはいったい幾つになるでしょうか？
［ 図1.4 ］砂地のある滑り台、その上からボールを転がす絵。
まず、このボールも質点として扱います。簡単の為、ボールを蹴っ飛ばした方向にｘ軸を取り、この軸だけに注目してみましょう。初期位置を原点とすると、時刻＜ｔ＞のときの距離は時間の関数ですので＜ｘ（t）＞となります。これに出題された条件を当てはめると、蹴っ飛ばした瞬間は＜ｘ（0）＝0＞、止まった瞬間は＜ｘ（5）＝10＞です。
さて、では速さを求めてみましょう。10メートルの距離に5秒掛かったわけですから、このボールは1秒間に2メートル進むペースで進んでいましたよ、ということが分かります。では、"このボールの転がっていたときの速さは、秒速2メートルです"と片づけてしまって良いのでしょうか？ボールの動きをよく観察すると、滑り台から砂地に落ちる直前のボールはかなり速く動いていましたし、逆に砂地では非常にゆっくり動いていました。砂地に落ちる直前のボールの速さのまま運動を続けたとしたら、1秒間に進む距離はもっと長そうな気がします。そうです、秒速2メートルという値はこのような瞬間々々の速さではなく、ボールの転がる運動全体の平均の速さでしかないわけですね。
瞬間の速さを調べるには、もっと計測間隔をもっと短くしてやれば良いのです。ボールが1秒の間に進んだ距離、0.1秒の間に進んだ距離を事細かく調べます。砂地に落ちる直前の0.001秒間に進んだ距離を測ってやり（どうやって測るかは永遠に秘密ですよ）、もし1センチメートル進んだと判明しました。その瞬間と同じ速さで1秒間進んだとすれば、0.001秒を1000倍してやり10メートルとなり、結果として秒速10メートルであるということが分かります。やはり、運動全体の平均である秒速2メートルより遙かに速いですね。ここから直感的にイメージできると思いますが、この0.001秒という時間間隔はさらに小さければ小さいほど、正確な瞬間の速さを見つけることができます。
［ 図1.5 ］ボールの位置を、時間を細かくした絵。
今の話を数式に表してみましょう。時間を時刻＜t＞から少し進めて＜t+Δt＞となったとき、ボールの位置は＜ｘ（t）＞から＜ｘ（t+Δt）＞まで動きました。この時、＜Δt＞だけ進んだ時間に対してボールの動いた位置差分＜Δｘ（t）＞は
［ 式1.1 ］Δｘ（t） = ｘ（t+Δt） - ｘ（t）
となり、＜t＞から＜t+Δt＞間の平均の速さは
［ 式1.2 ］Δｘ（t） / Δt
と表すことが出来ます。このときの＜Δt＞を限りなく0に近づけるために、極限記号＜lim［x→a］＞を用いて
［ 式1.3 ］ｖx（t） = lim［t→0］（ Δｘ（t） / Δt ）
と書くことで、これを時刻＜t＞における瞬間の速さと呼ぶことが出来ます。
このΔを、限りなく0に近いという意味で＜d＞という記号で表すと以下のようになります。
［ 式1.4 ］ｖx（t） = dｘ（t） / dt
このような＜dｘ（t） / dt＞、つまり＜ｖx（t）＞を、＜t＞に関する＜ｘ（t）＞の導関数と呼びます。そして＜ｖx（t）＞を求めることを、＜dｘ（t）＞を＜t＞で微分（一階微分と呼びます）するといいます。もっと具体的な言い方をすれば、位置＜ｘ（t）＞を時刻＜t＞で微分した値が、速さ＜ｖx（t）＞なのです。
式1.4からイメージできるように、＜ｖx（t）＞は＜ｘ（t）＞が微少時間の間にどれだけの変化が生まれるのか、を表しています。これを"時間変化率"といいます。
以上の話は位置ベクトルの1要素のみで話を進めてきましたが、これを全ての要素に当てはめてみると、
［ 式1.5 ］
ｖx（t） = dｘ（t） / dt
ｖy（t） = dｙ（t） / dt
ｖz（t） = dｚ（t） / dt
となります。これらをまとめて
［ 式1.6 ］ｖ（t） = dｒ（t） / dt
となり、この＜ｖ（t）＞を速度ベクトルと定義します。速度ベクトルは位置ベクトルと違って、座標系の基準点の取り方によって値が変わるものではありません。
なお、細かいことですが、物理世界では"速さ"と"速度"を明確に別のものであるとしています。速さといった場合、それは大きさだけで方向を持たないスカラー量を指します。速度といった場合、それは大きさと方向を持ったベクトル量を指します。例えば同じ速さで円運動をしている質点があった場合、その速度は一定ではありません。なぜならば、この質点の向かっている方向は円の接線方向に絶えず変化しているからです。よって速度も時間によって絶えず変化しているのです。
［ 図1.6 ］円運動をしている質点。その接線方向に速度ベクトルが引かれている絵。
1.1.4 加速度.
＜ｖ（t）＞という書き方からも分かるとおり、速度も位置と同じように"時間によって変化する"時間の関数です。位置が時間によってどう変わるかを表記するために速度を定義したように、速度が時間によってどう変わるかを見ていきます。前項で説明した道を同じように辿ってみますと、時刻＜t＞における速度の時間変化率＜ａ（t）＞は以下のようになります。
［ 式1.7 ］ａ（t） = lim［t→0］（ Δｖ（t） / Δt ）
= dｖ（t） / dt
この速度の時間変化率＜ａ（t）＞を"加速度"と呼びます。私たちは、電車が走っている最中よりも発進/停止時に違和感を覚えます。このように、肌で感じることの出来るものはむしろ速度よりもこの加速度であるといえます。
式1.7は式1.6と組み合わせて、
［ 式1.8 ］ａ（t） = d^2ｒ（t） / dt^2
と書くこともできます。この式から分かるとおり、加速度は位置を２回だけ微分（二階微分と呼びます）した時間の関数です。では、この加速度、すなわち質点の運動を変化させる要因はどこから生まれるのでしょうか？この部分について、1.2節にて詳しく説明していきます。
1.1.5 積分による変位計測.
前項までの話によって、質点の位置の時間による変化から、速度及び加速度を求めることが出来ました。今度は逆に、それぞれの時刻での速度や加速度が分かっている場合に、質点の位置がどう変化するのかをみていきます。このような、位置の変化分を変位と呼びます。
［ 図1.7 ］発進から１秒ごとの速度が表に並んでるグラフ.
まずは図1.7のような、車が発進するときの速度表があったとします。発進から最初の0.5秒間目は、秒速0.01メートルで進んでいたと図から読みとれます。ここから、0秒目から1秒目に進んだ距離はおおよそ0.01メートルじゃないかと推測できます。同様に1.5秒目が秒速0.05メートルですので、1秒目から2秒目に進んだ距離はおおよそ0.05メートル、と計算していきます。
このように一定時間毎にどの程度の距離を進んだのか記録しておいて、最後に全ての距離を足してやれば、それが車の進んだ距離としてそれなりに正しい値が出せるでしょう。これは図1.8のように、単位時間毎の速度と微少時間を掛けたものを面積として、足し合わせたものとイメージすることが出来ます。
［ 図1.8 ］図1.7を１秒区切りで、その部分を長方形に切り取って加算！な絵。
図1.8では１秒刻みで計算されていますが、この時間の間隔をもっと細かくしてやることで、さらに正確な距離が求まるであろうことは想像に難くありません。実はこのことは、式1.6の両辺に＜dt＞を掛けてやることで表現可能です。
［ 式1.9 ］ｖ（t） * dt = dｒ（t）
式1.9は、"速度に微少時間をかけたもの"が"微少時間における位置の変化量"であることを意味しています。この変化量を積み重ねていくことで、積み重ねた時間分の位置の変化量が求まります。これが微分の逆演算、"積分"の基本的な考え方です。
微分のときと同じように、時間の間隔を微少時間＜Δt＞として、時刻＜t＞での瞬間の速度を＜ｖ（t）＞とすると、時刻＜t1＞から時刻＜t2＞までに車の進んだ距離＜ｓ＞は
［ 式1.10 ］ｓ = ｖ（t1）Δt + ｖ（t1+Δt）Δt + ... + ｖ（t2）Δt
= Σｖ（ti）Δt
と表記されます。（Σ記号は総和を表します）このとき、＜Δt＞を極限まで０に近づけた
［ 式1.11 ］ｓ = lim［Δt→0］（ Σｖ（ti）Δt ）
こそが、正しい距離を表しているといえます。これでは式が見づらいですので、積分記号＜∫＞と、微分でも用いた＜dt＞を用いて
［ 式1.12 ］ｓ = ∫［t1-t2］ｖ（t）dt
と書くことにしましょう。この書き方で＜t1＞から＜t2＞までの全ての＜ｖ（t）dt＞を加算したことになり、この間に動いた変位を出すことができます。時刻＜t1＞での位置が予め分かっていれば、時刻＜t2＞での位置も求まることになります。これが積分と呼ばれる計算です。微分が瞬間の速度を求めているのに対して、積分はその瞬間の速度を一定時間ひたすらかき集めているのです。
以上は位置と速度の関係について話を進めていましたが、これは速度と加速度の関係であっても同様の式が立てられます。この積分を使うことで、加速度から速度が求まり、速度から位置が求まるのです。
1.2 質点の運動力学.
1.2.1 ニュートンの三運動法則.
質点の運動には、よく知られた三つの"ニュートンの運動法則"が成り立つことが分かっています（※）。古典力学において全ての運動は、この三つの運動法則から導かれることになります。どれもよく考えれば直感的に理解できる法則ですので、覚えるより理解することを優先して読み進めてみてください。
（※） 非慣性系では成立しないが、ここではその話には触れないことにする。
1.2.2 慣性の法則.
［ニュートンの運動第一法則：慣性の法則］
他からの影響が無いとき、質点は現在の運動を継続する。
［ 図1.9 ］ボールを転がすだけの絵。
以下のような実験を考えてみましょう。図1.9のようにボールを平らな地面に転がしてみると、Ａ地点でボールの動きは止まりました。地面をぴかぴかに磨いて摩擦を減らした後で、もう一度同じように転がしてみました。すると、今度はＡ地点よりも遠いＢ地点で止まりました。さらに摩擦を減らすと、さらに遠いＣ地点で止まりました。この実験から推測されるのは、もし摩擦という影響がボールに働かなければ、ボールは最終的に地の果てまで転がり続けるのではないか？言い換えれば、何かしらの影響が働かない限り、質点は同じ動きをし続けるのではないか？ということです。これを法則として定めたのが、ニュートンの運動第一法則、よく聞く言い方で"慣性の法則"です。この法則は、今後の話を進める上での大前提となります。質点は、他からの影響が無い限り、静止していれば静止を続けますし、運動していれば等速直線運動を続けるのです。
運動を継続する、というのは、1.1節で説明したとおり「加速度は０のままである＝速度が変わらない」と言い換えることができます。では、"他からの影響"とはいったい何物なのでしょうか？……おおよそ察しがついていますね、そう、"力"です。"力"とは要するに、"質点の運動を変化させる何かしらの影響"を、物理用語でそう呼称しているものなのです。先ほどの話で出た"摩擦"や、"物体同士の衝突時の跳ね返り"といった物体の衝突・接触時に起こる運動の変化も、"重力"といった非接触時の運動の変化も、全部引っくるめて"力が働いた"と呼称します。
ボールを全力で上空に投げ上げても、いつかは地面に落っこちてきてしまいます。これは重力や空気抵抗といった力がボールに対して働いているからであって、もしもこれらの力が働いていなければ、ボールはそのまま宇宙の果てまで等速直線運動で飛んでいってしまうことでしょう。
1.2.3 慣性質量.
［ 図1.10 ］段ボール箱二つ。片方はぱんぱんに荷物が入っている図.
図1.10のように、荷物がたくさん入った、重たい段ボール箱と、中身が空の状態の、軽い段ボール箱を、両方とも持ち上げてみるとします。このとき、重たい段ボール箱が動かす方がより体力がいるのは自明ですね。段ボール箱の運動を変化させるときに、重たいか軽いかでどれだけの体力を必要とするかが違ってしまっているのです。つまり、物体にはそれぞれ"運動の変化に抵抗するパラメータ"が存在するのです。私たちは普段このパラメータを"重量"と呼んでいます。
ここで一つ問題があります。困ったことに、私たち人間は遠い未来、月面や火星といった遠い星々で暮らすことになるでしょう。しかし、月面では地球よりも重力が弱く、同じ重さの物体を持ち上げても地球ほどは体力を要しないでしょう。つまり"重量"は環境によって変わってしまうのです。では重力の全くない宇宙空間ならば、どんなに重たい物体、例えばスペースシャトルであっても人間の手で軽々と動かせるかというと、勿論そんなことはありません。重力に逆らって持ち上げた時に重力の分だけ軽く感じるのであって、左右（左右という呼び方が適切かはともかく）に物体を動かすときは、地球上でも月面上でも宇宙空間でも同じだけの体力が必要なのです。例えば図1.10がもし宇宙空間で重力が全くなかったとしても、非力な小指だけを使って重たい方の段ボールは持ち上げられないでしょう。
このような、重力に左右されない"物体固有の、大きいほど運動を変化させづらくなる量"を"慣性質量"もしくは"質量"（※）といい、通常は＜ｍ＞という記号で表します。そして重量と明確に区別するために「大きい」、「小さい」といった言い方をします。質量と重量は比例の関係にあります。地球上では質量が大きい物体は重たいですし、質量が小さい物体は軽いのです。
（※） 厳密には質量には"慣性質量"と"重力質量"があるが、これはどちらも同じ値であることが証明されている。
この質量を用いることで、将来地球を離れる未来の人類に対しても安心して今後の話を進めることができます。なお古典力学では、質量は時間によって変化することのないスカラー（実数）値として扱います（※）。
（※） 例えば相対論を考慮に入れると、質量は光速に近づくに連れて大きくなる。結果、慣性が強く働くようになり、最終的にはどれだけ加速しようとしても慣性のせいで速度が上がらなくなる。物体の速度が光速を越える事は出来ない、と言われているのはこのせいである。とはいえ、余程光速に近いレベルの速度を扱うのでなければ誤差の範疇に収まるので、事実上"質量は不変"として扱っても問題は無い。
1.2.4 ニュートンの運動方程式.
［ニュートンの運動第二法則：運動方程式］
質点に働く力は、質点の運動量の時間変化率に等しい。
では"力"を、もう少しきちんと整理してみましょう。その為にまずは"運動量"という物理量について考えます。運動量＜ｐ＞とは、質点の質量＜ｍ＞と速度＜ｖ＞を掛けたもので、質点がどれだけ運動に勢いがあるのかを表したものです。
［ 式1.13 ］ｐ = ｍｖ
ここで速度＜ｖ＞は時間の関数であり、実際には＜ｖ（t）＞であること、運動量＜ｐ＞も中身は速度にスカラー値の質量＜ｍ＞を掛けただけの時間の関数＜ｐ（t）＞であることを忘れずにおいてください。
同じ速度で運動している物体があったとき、質量が小さい物体よりも大きい物体の方が大きい運動をしている、といっても違和感はありませんね。同様に、質量が同じ物体があったとき、速度の遅い物体よりも速い物体の方が大きい運動をしている、というのも違和感ありませんね。そこで、質量と速度を掛けた値を運動量と呼び、運動量が大きければ大きいほど凄い運動をしている、としたわけです。1kgのボールを秒速4mで投げたときと、2kgのボールを秒速2mで投げたときの、両方のボールの持っている運動量は同じであるというのです。
［ 図1.11 ］ボール二つが飛んでいる。同じ運動量。の図。
ではなぜ、運動量などという単位を作ったかというと、力とは「運動量を単位時間当たりにどれだけ変化させることができるか」を表したものだからなのです。つまり、力＜Ｆ＞を運動量の時間変化率として
［ 式1.14 ］dｐ/dt = dｍｖ/dt = Ｆ
と定義するわけです。ここで、質量＜ｍ＞は時間による影響を受けませんから、
［ 式1.15 ］ｍdｖ/dt = Ｆ
と外に追い出すことができます。これは速度＜ｖ＞で表現されていますが、以下の式のように加速度＜ａ＞や位置＜ｒ＞で表現しても意味は変わりません。
［ 式1.16 ］ｍａ = Ｆ
［ 式1.17 ］ｍd^2ｒ/dt^2 = Ｆ
これが、目にすることも多いでしょう**"ニュートンの運動方程式"**です。
慣性の法則のとおり、力が働いていなければ加速度が０であるということが確認できます。質量＜ｍ＞はスカラー値なので、ベクトル値である＜Ｆ＞と＜ａ＞は同じ方向を向いている、ということを式1.16は暗に示しています。つまり物体の加速は、力の働いた方向にそのまま生まれるのです。質点の項で確認したとおり、質量＜ｍ＞が大きければ大きいほど、同じ力＜Ｆ＞でも加速度＜ａ＞が小さくなる、ということもこの式では表されています。
1.2.5 微分方程式.
式1.16によって、任意の時刻＜t＞に質点に働く力＜Ｆ＞（※）が分かれば加速度＜ａ（t）＞を求めることが出来るようになりました。
（※） この力＜Ｆ＞は時間や質点の位置など、様々な要因が絡まった関数であると考えられる。ので正式には＜Ｆ（t,ｒ（t）,...）＞となる。
しかし、質点の運動、つまり位置や速度はこの段階ではまだ分かっていません。式1.17のような、未知な関数（ここでは位置）の微分形を含む方程式を、一般に"微分方程式"と呼びます（※）。
（※） 二階微分で、変数が時刻＜t＞一つであることから、正確には"二階の常微分方程式"である。
ゲームの場合、この微分方程式に対して初期条件は揃っていることと思います。初期条件とは、時刻＜t＞が0の時の位置ベクトル＜ｒ（0）＞や速度＜ｖ（0）＞のことを指します。この初期条件と式1.17のような微分方程式から、任意の時刻＜t＞における位置ベクトル＜ｒ（t）＞や速度＜ｖ（t）＞を求めることを初期値問題といい、本章での最終目標となります。実際の解法については、1.3節にて詳しく説明します。
1.2.6 作用反作用の法則.
［ニュートンの運動第三法則：作用反作用の法則］
二つの質点に対して、それぞれに働く力の和はゼロになる。
この法則は、前述された二つの法則とは異なり、二つの質点が与える関係を記したものになります。例えば机に置かれた本は、本が重力によって落ちようとしています、が、机はその力を押し返しています。本と机がそれぞれに対して働いている力は、大きさが同じで正反対の方向を向いています。
［ 図1.12 ］机の上に置かれている本。その力関係。
本という物体の中で、重力によって発生している力と、机から押し返された力が釣り合い、結果として運動量は変えるには至らず、位置はそのまま机に"置かれている"わけです。この机と本をそれぞれ質点と考え、関係を数式で表すと、質点iがjに与える力を＜Ｆij＞としたとき
［ 式1.18 ］Ｆij = - Ｆji
となります。このとき力はベクトルなので、力の大きさが等しいというだけではなく力の働いている向きが正反対であることも表しています。この法則は、第二章の"質点系"を学んだ時点で、"運動量保存則"という法則に言葉を変えて、いかんなく効果を発揮します。
1.2.7 角運動量.
角運動量とは、ある基準点から見た質点の回転に対する勢いを定義したものです。質点の基準点からの位置ベクトルを＜ｒ＞、運動量を＜ｐ＞としたとき、質点の（基準点からみた）角運動量＜ｌ＞は
［ 式1.19 ］ｌ = ｒＸｐ
となります。
［ 図1.13 ］角運動量を指図した絵1。
運動量は座標系によって大きさは変わりませんが、角運動量は座標系の原点の取り方に依存することに注意してください。てこの原理からイメージできるように、原点から離れているほど、運動量が大きいほど角運動量は大きくなります。さらに、位置ベクトルの向き方向への速度は回転成分と見なされないので、速度の方向が位置ベクトルの向きに対して垂直であるほど角運動量は大きくなります。
角運動量は運動量を原点周りの回転の勢いとして変形させただけのものであり、運動量とは全然別とする"物体の回転を表したパラメータ"ではないことに注意してください。
［ 図1.14 ］角運動量を指図した絵2。
ここで、式1.19の両辺を時間で一階微分すると以下のようになります。
［ 式1.20 ］dｌ/dt = d（ｒＸｐ）/dt
= dｒ/dt×ｐ + ｒ×dｐ/dt
= ｖ×ｐ + ｒ×Ｆ
= ｖ×ｍｖ + ｒ×Ｆ
＜ｖ×ｍｖ＞は外積の性質より０になりますので、
［ 式1.21 ］dｌ/dt = ｒ×Ｆ
となります。この＜ｌ＞の時間変化率＜dｌ/dt＞を、トルク＜ｔ＞（※）と呼びます。
（※） トルクは力のモーメントとも呼ばれる。
［ 式1.22 ］ｔ = ｒ×Ｆ
この式から分かるとおり、トルクとは、力が働いたとき基準点から見て回転方向にどれだけの勢いがつくかを表しています。この角運動量及びトルクは第二章にて大いに活躍します。
1.3 微分方程式の初期値問題.
1.3.1 微分方程式の初期値問題.
1.1節にて、物体の運動に対して速度と加速度がどう定義されているのかを説明しました。1.2節では、物体に働く力と、それによって発生する加速度の関係を導きました。後はこの加速度から、質点がどのように位置を変化させていくのかを記述することが出来れば、今までの話が全て繋がることになります。その為には、項1.2.5で説明したとおり、式1.17の微分方程式を解いてやる必要があるわけです。
解法は、大別して"解析解法"と"数値解法"の二つに分けられます。このうち解析解法は条件が非常に厳しく、ゲームには殆どのケースに於いて使用できない、と考えてしまって問題ありません。そこで数値解法に頼ることになるのですが...その辺りを含めて、解説をしていきます。
1.3.2 解析解法.
説明が煩雑になりそうなので、実際に具体例を出して話を進めていきます。
［ 図1.15 ］放物運動の絵.
地面からの高さ＜ｈ＞の場所からx軸方向へ、地表からの角度＜θ＞、速さ＜ｖ0＞で重さ（質量）＜ｍ＞のボールを投げた場合の運動について考えていきます。ボールには一定の重力が働いているとします。これは放物運動と呼ばれ、質点の運動の中でも最もシンプルでよく知られたものです。
初期条件の成分表示は問題から、以下のようになります。
［ 式1.23 ］ｒx（0） = 0.0, ｒy（0） = ｈ, ｒz（0） = 0.0
ｖx（0） = ｖ0cosθ, ｖy（0） = ｖ0sinθ, ｖz（0） = 0.0
重力については三章にて説明を行いますが、今のところは重力加速度を＜ｇ＞として
［ 式1.24 ］Ｆ = ｍｇ
が鉛直方向へ働く、としておきます。ここから、ボールの運動方程式の成分表示は以下のように立てることが出来ます。
［ 式1.25 ］ｍd^2ｒx/dt^2 = 0
ｍd^2ｒy/dt^2 = - ｍｇ
ｍd^2ｒz/dt^2 = 0
これらの式を積分し、初期条件を利用すると
［ 式1.26 ］dｒ（t）x/dt = ｖ0cosθ
dｒ（t）y/dt = - ｇt + ｖ0sinθ
dｒ（t）z/dt = 0.0
となり、この式がボールの速度を表しています。さらに積分すると
［ 式1.27 ］ｒ（t）x = ｖ0cosθt
ｒ（t）y = - 1/2ｇt^2 + ｖ0sinθt + ｈ
ｒ（t）z = 0.0
と、ボールの位置が式で表せました。この式を用いることで、任意の時刻におけるボールの位置を求めることができます。ボールが落ちる地点が原点からどれだけ離れているかを求めることも、どの角度＜θ＞で投げるのが一番遠くに飛ぶのかも、式1.27があればできますよね？このように、積分を用いて解析的に微分方程式を解く解法を解析解法と呼びます。
解析解法は、計算の精度によらず、好きな時刻＜t＞に対する答えが（式に＜t＞の値を代入するだけで）必ず見つかるのが大きな利点です。しかし、解析解法には致命的な欠点があります。それは、式を立てて解ける運動が非常に限られていることです。
今度は先ほどの放物運動の話に、ボールに対して空気抵抗が働く場合のことも考えてみましょう。初期条件は変わりません。空気抵抗については重力と同様、三章にて説明を行います。が、今のところは粘性抵抗定数を＜ｋ＞、ボールの速度を＜ｖ＞として
［ 式1.28 ］Ｆ = ｋｖ（t）
が速度の反対方向に働く、としておきます。つまり、ボールが速ければ速いほど、それに反発する力が働くわけです。先ほどと同じようにボールの運動方程式を立ててみると、
［ 式1.29 ］ｍd^2ｒx/dt^2 = - ｋｖx（t）
ｍd^2ｒy/dt^2 = - ｍｇ - ｋｖy（t）
ｍd^2ｒz/dt^2 = 0
となります。この式の積分を求めるのは不可能ではありませんが、そう簡単にはいきません。「これくらいなら解けるよ」という方も、もし"ボールを投げ上げて10秒後に、カラスにぶつかって明後日の方向に落ちていく"という条件が付いたとしたら、式を立てろと言われてもお手上げでしょう。少し複雑な条件が付いただけでも、解析解法を利用するのは非常に厳しくなってしまうのです。
1.3.3 数値解法.
このような場面で役に立つのが数値解法です。解き方は以下のようになります。
1.1.5項で説明した積分の計算を思い出してみてください。微少時間毎の速度をかき集めて、少しずつ位置の変位を更新していく計算を行っていましたね。この微少時間を積分計算では無限小時間であるとしていますが、これを有限時間にしたとしても、それなりに正しい答えは出てきそうです。つまり、積分では連続時間（無限に小さな時間が流れている世界）を元に計算を行いますが、数値解法ではこれを離散時間（有限に細かい時間が流れている世界）で扱うわけです。当然、この有限時間は小さければ小さいほど、本来の積分計算より誤差は少なくなります。
例えば、時刻＜t＞での質点の速度が＜ｖ＞で、位置が＜ｒ＞であったとします。ここから有限時間＜Δt＞進めて＜t+Δt＞としたとき、速度と位置はどうなっているでしょうか？これを解くことができたとしたら、そこから少し時間を進めた＜t+Δt+Δt＞も、さらに進めた＜t+Δt+Δt+Δt＞も同じように解くことができそうです。このように、少しずつ時間を進めながら質点の運動を追いかけていくやり方が、数値解法と呼ばれる解法です。
この考え方を元にして、近似式を立てると以下のようになります。
［ 式1.30 ］ａ（t） ≒ { ｖ（t） - ｖ（t-Δt） } / Δt
ｖ（t） ≒ { ｒ（t） - ｒ（t-Δt） } / Δt
これはあくまで近似式であり、他の近似式を立てることも可能です。
この為、数値解法には様々なアルゴリズムが生まれます。
［ 式1.31 ］ａ（t） ≒ { ｖ（t+Δt） - ｖ（t） } / Δt
ｖ（t） ≒ { ｒ（t+Δt） - ｒ（t） } / Δt
ここから、微分方程式を近似した方程式（差分方程式といいます）を立てると
［ 式1.32 ］ｖ（t+Δt） ≒ ｖ（t） + ａ（t）Δt
ｒ（t+Δt） ≒ ｒ（t） + ｖ（t）Δt
となります。この差分方程式の立て方にはいくつかの方法があり、式1.32の方法はオイラー法と呼ばれる最もポピュラーなものです。この差分方程式を立てることが出来れば、後は式1.29と組み合わせることによって、初期時刻からの差分を累算していくことによって任意時刻での位置や速度を求めることができます。
数値解法は、解析解法と違い式を解いたりする必要が無く、ゲームで用いるのには適した解法といえます。実際、式1.32をオイラー法と知らず知らずのうちにゲーム中のキャラクタ制御に使っていたという方は多いでしょう。
ただし、数値解法はあくまで近似式であり、有限時間＜Δt＞の取り方で結果が変わってしまう、ということを忘れてはいけません。＜Δt＞が小さければその分処理速度を喰われることになり、大きければ計算誤差として跳ね返ってきます。最悪なのはこの誤差が指数関数的に跳ね上がり無限大に近づく、俗に"発散する"と呼ばれる状況になってしまうことです。数値解法を用いる限り、この発散に関してはずっと頭を悩ませる要因となります。
chapter 2 剛体の運動力学.
前章では、質点と呼ばれる大きさを無視した物体のことを考え、これを用いて物体の運動について考えていきました。打ち上げられた球やパーティクルの計算など、形状を考慮にいれなくても構わない、位置のみが意味を持つ状況下では、この質点の運動力学を用いることで一連の運動を掌握しプログラミングすることが出来ました。
本章では、大きさを持った物体を取り扱います。ただし、いきなり"流れる水"や"人間"といった複雑なものに取りかかれるわけではありません。まずは物体の中でも一番単純な部類である"剛体"について話を進めていきます。
［ 図2.1 ］木箱が積まれている絵.
剛体とは、図2.1の木箱のような、大きさを持ち変形することがない固体のことをいいます。大きさがあると言うことは、形があり、向きがあり、回転します。質点と唯一違う点はここです。そして壊れたり、凹んだりすることがありません。常に同じ形状をしている、という制限を設けることによって、計算を簡単にしています。それ故、剛体は、形状があって硬いオブジェクトの挙動計算などに適しています。オブジェクトが壊れたり変形したりといった高望みをしなければ、ゲーム中で使用するオブジェクトの殆どはこのタイプに属するのではないでしょうか。
質点と剛体は"大きさがある"というだけの違いしかありませんが、それでも前章と比べても難解な式が並ぶことになります。しかし、実際問題としてやっていることはニュートンの運動方程式をこねくり回しているだけであって、それほど高度なことをしているわけではありません。まずは剛体を理解するための土台として、"質点系"という概念について解説を行い、その後で剛体について話を進めていきます。
2.1 質点系.
2.1.1 質点系とは.
二個以上の質点を一括りにして物事を考えるとき、これを質点系と呼びます。括り方は観測者、つまりあなたの自由です。例えばボーリングのシミュレーションを計算しようとして、ピンを質点としたならば、ボーリング場の１レーンに並べられたピンをまとめて質点系として考えるのがスマートでしょう。隣のレーンのピンを含めて質点系としてしまっても問題はありません...が、隣のレーンの状況が目の前のピンに対して影響を与えるとは考えにくいので、この場合は無視してしまった方が都合が良いですよね。計算しやすい、分かりやすい括り方を行うのが良いでしょう。
［ 図2.2 ］複数個のボールが袋の中に入ってる図。それぞれに位置ベクトルが存在する.
まず、何個かのボールを袋の中に入れます。このボールを質点として考えて、袋の中のボールの集まりを質点系とします。この袋の口を縛って、思いきり力を込めて上空に放り投げました。このとき、中のボールどれか一つに着目して動きを観察してみると、他のボールにぶつかったり袋に引っ張られたり、滅茶苦茶な軌跡を描いているのが分かります。他のボールも同様に、規則性のある動きはしていないのは想像に難くありません。どのボールも非常に複雑な運動で、これを式に立てろと言われても頭を抱えてしまいます。
しかし、袋全体を一つの質点として見てみると、特に複雑な動きをしているわけでもなく、単純な放物線に沿って動いているようです。袋自体の形はグニャグニャと変わっているでしょうけれども、今回は袋全体が質点であり、袋の中心位置だけに注目しています。不思議なことに、この袋のどこかにある中心位置の軌跡を辿ると、一つのボールを放り投げた時と同様に放物運動をしているのです。
［ 図2.3 ］袋を投げたら放物運動だった！絵。
ここで一つの疑問点が生まれるでしょう、"袋の中心位置とは一体どこを指しているの？"と。これは袋の中のどのボールの中心でもなければ、袋に適当に目印を書き込んでおいた場所でも、締めた袋の口でもなさそうです。そもそも、ボールを質点と見立てたときの中心位置、という場所ですら非常に曖昧です。前章ではボールの中心位置に点を取ってそこを質点として運動を見ていきましたが、なぜ中心位置が使われたのでしょうか？ボールならまだしも、ボーリングのピンを質点としたときの中心位置はどこにしましょうか？と聞かれて即座に答えられる人は少ないでしょう。この"物体の中心位置"とは一体なんなのでしょうか？
2.1.2 運動量保存則.
中心位置を探る前に、まずは袋の中のボールについて考えていきます。それぞれのボールを質点として、各質点＜i＞は質量＜ｍi＞と位置ベクトル＜ｒi＞で指し示すとします。これらのボールをまとめて質点系として括ります。質点系内部の質点同士で働く力を"内力"といい、質点＜i＞が質点＜j＞に受ける内力を＜ｆji＞と記します。それ以外の、系の外から各質点が受ける力を"外力"といい、質点＜i＞の受ける外力を＜ｆi＞と記します。今回のような袋の放り投げケースでは、外力には重力や空気抵抗などが挙げられます。同様に内力は、ボール同士の衝突や摩擦などが挙げられます。
［ 図2.4 ］袋の中のボール二つ。重力を矢印で書いて外力として書き、衝突の力を矢印で書いて内力と記す。
袋の中のボールの個数をN個として、各質点の運動方程式（式1.16）を書き記すと以下のようになります。
［ 式2.1 ］ｍ1ａ1 = ｆ1 + ｆ21 + ｆ31 + ... + ｆN1
ｍ2ａ2 = ｆ2 + ｆ12 + ｆ32 + ... + ｆN2
...
ｍNａN = ｆN + ｆ1N + ｆ2N + ... + ｆ（N-1）N
これらの式の左右の辺を足し合わせてみると、内力＜ｆij＞には対応する＜ｆji＞が必ず存在することが分かります。この二つの内力は前章1.2.6項で習ったとおり、作用反作用の法則が適用されます。
［ 式2.2 ］ｆij = - ｆji
よって全ての内力＜ｆij,ｆji＞は打ち消しあい、結果として
［ 式2.3 ］ｍ1ａ1 + ｍ2ａ2 + ... ｍNａN = ｆ1 + ｆ2 + ... ｆN
Σｍiａi = Σｆi
という式に落ち着きます。このうち、右辺は各質点に働く外力の総和を表しており、これを＜Ｆ＞と置き換えます。
［ 式2.4 ］Σｍiａi = Ｆ
質点＜i＞の運動量＜ｐi＞は、速度＜ｖi＞として
［ 式2.5 ］ｐi = ｍiｖi
であり、全質点の運動量の総和を＜Ｐ＞として
［ 式2.6 ］Ｐ = Σｐi
と置き換えるとすると、式2.4から式2.6を組み合わせて以下の式に辿り着きます。
［ 式2.7 ］Ｆ = Σｍiａi
= Σｍidｖi/dt
= dΣｍiｖi/dt
= dΣｐi/dt
= dＰ/dt
式2.7を直訳すると、"全質点の運動量の一階微分、つまり時間変化率は、全外力の総和に等しい"となります。右辺の値が０、"質点系に対して外力が一切働いていない状態"ならば、質点系全体の運動量は一定を保たれるのです。質点系内部で衝突や摩擦などの内力が発生して、それぞれの質点の運動状態は変化したとしても、その運動量の合計は変わらないままなのです。このような法則を"運動量保存則"といいます。
式がずらずらと並んでばかりで鬱陶しいので、ペットボトルに空気をギュウギュウに詰めて飛ばす、俗に言う"ペットボトルロケット"を例を挙げて運動量保存則を説明しましょう。このペットボトルロケット、実は空気だけよりも、水を少し入れておいた方が高く飛ぶのですが、この理由を運動量保存則で簡単に説明ができます。
［ 図2.5 ］ペットボトルロケット二つ.片方は水入り.
まず、ペットボトルと、ペットボトル内の空気と水を質点系として括ります。蓋をして静止している状態では、それぞれの運動量は０のまま一定の状態を保っています。ここで蓋を外してみると、圧縮された空気の圧力によってペットボトルの口から空気が逃げようとします。この時空気は、作用反作用の法則によりペットボトル本体に対して逆方向に力を働かせます。これが推進力となり、ペットボトルは宙に飛ぶわけです。しかし、空気の質量はペットボトルよりかなり小さいので、空気の運動量は対した大きさではなく、その反作用である推進力もそれほど大きくなりません。この時もしペットボトルに水が入っていたらどうなるでしょうか。空気の進行方向に、空気より遙かに大きい質量の物体が飛び出すことになります。つまり運動量は非常に大きなものとなり、ペットボトルの受ける推進力も非常に大きくなるのです。故に水を入れたペットボトルの方がより高く飛ぶわけです。
［ 図2.6 ］ペットボトルロケット発射絵.力と運動量のベクトルを記す.
2.1.3 質量中心.
それでは、いよいよ物体の中心位置について触れることにしましょう。質点系での全質点の合計質量を＜Ｍ＞とします。先ほどの話で言えば、ボール全部の質量と袋の質量を足したものが合計質量です。
［ 式2.8 ］Ｍ = ｍ1 + ｍ2 + ... + ｍN
= Σｍi
ここで、質量中心（※）＜Ｒ＞を式2.9のように定義し、その一階微分、二階微分がそれぞれ式2.10、式2.11であるとします。（※） 重心とも呼ばれるが、本文章では質量中心で用語を統一する。
［ 式2.9 ］Ｒ = Σｍiｒi ／ Ｍ
［ 式2.10 ］dＲ/dt = Σｍiｖi ／ Ｍ
［ 式2.11 ］d^2Ｒ/dt^2 = Σｍiａi ／ Ｍ
この質量中心＜Ｒ＞を利用すると、質点系全体の運動方程式を表していた式2.7は
［ 式2.12 ］Ｆ = Ｍd^2Ｒ/dt^2
と書き換えることが出来ます。
この式をよくよく見ると、質点の運動方程式である式1.17が"位置ベクトルｒ＝質量中心Ｒ"、"質量ｍ＝合計質量Ｍ"として置き換えられただけであることが分かります。つまり、複数の質点として扱っていたはずの質点系に対して質量中心＜Ｒ＞と合計質量＜Ｍ＞を求めてやることで、まるで一つの質点を扱うかのように運動を扱うことができるのです。一章で"物体を点と考えて質点とみなし、運動を観測する"ことが許されたのは、暗黙的に物体の質量中心＜Ｒ＞と合計質量＜Ｍ＞を扱っていたからに他ならないのです。
式2.12を用いれば、物体がどんなに複雑な形状をしていた場合であっても、その質量中心と合計質量が求まれば後は１つの質点として扱うことができてしまいます。先ほどのボールが入った袋でいえば、袋全体に対する質量中心である位置が放物運動の軌跡を描いていたのです。
［ 図2.7 ］袋には質量中心と合計質量があった！それが放物運動の軌跡を描いていたんだ！絵.
2.1.5 角運動量保存則.
質点系での各質点の角運動量＜ｌi＞とトルク＜τi＞も注目してみましょう。
［ 式2.13 ］ｌi = ｒi Ｘ ｐi
［ 式2.14 ］τi = ｒi Ｘ ｆi
一章で既に解説しましたが、再度確認をしておきます。角運動量とは、ある基準点から見た質点の回転に対する勢いを定義したものです。そして角運動量の時間変化率（一階微分）がトルクです。運動量に比べて今ひとつイメージし辛いかもしれませんが、今後の話の流れにおいて非常に重要な意味を持つ値ですので、しっかりと定義を頭にたたき込んでおいてください。
全質点の角運動量の総和を＜Ｌ＞とすると、以下のようになります。
［ 式2.15 ］Ｌ = Σ（ ｒi Ｘ ｐi ）
各質点に働く外力と内力を式2.14の右辺＜ｆi＞に入れると、以下のような式になります。
［ 式2.16 ］τi = ｒi Ｘ ｆi + ｒi Ｘ ｆji + ｒi Ｘ ｆki + ... + ｒi Ｘ ｆNi
これを各質点毎に立ててみると、
［ 式2.17 ］τ1 = ｒ1 Ｘ ｆ1 + ｒ1 Ｘ ｆ21 + ｒ1 Ｘ ｆ31 + ... + ｒ1 Ｘ ｆN1
τ2 = ｒ2 Ｘ ｆ2 + ｒ2 Ｘ ｆ12 + ｒ2 Ｘ ｆ32 + ... + ｒ2 Ｘ ｆN2
...
τN = ｒN Ｘ ｆN + ｒN Ｘ ｆ1N + ｒN Ｘ ｆ2N + ... + ｒN Ｘ ｆ（N-1）N
となります。ここで、2.1.2項の時と同じように左右の辺を足し合わせます。このとき、作用反作用の法則（式2.2）によって
［ 式2.18 ］ｒi Ｘ ｆji + ｒj Ｘ ｆij =
ｒi Ｘ ｆji - ｒj Ｘ ｆji =
（ｒi - ｒj） Ｘ ｆji
と変形でき、＜（ｒi - ｒj）＞と＜ｆji＞が同じ方向を向いていることから、この式の外積の結果は0になります。従って、＜ｒi Ｘ ｆji＞と＜ｒj Ｘ ｆij＞は打ち消しあって消えてしまいます。その結果、外力によるトルクだけが式の中に残ることになります。全質点のトルクの総和を＜Τ＞と置き換えるとすると、
［ 式2.19 ］dＬ/dt = Τ = Στi = Σ（ ｒi Ｘ ｆi ）
となります。この式2.19から分かるとおり、質点系での全質点のトルク（角運動量の時間変化率）は、全外力によるトルクの総和に等しいということになります。つまり、最右辺の値が０、"質点系に対して外力のトルクが一切働いていない状態"ならば、質点系全体の角運動量も一定に保たれるのです。質点系内部でどんな内力によるトルクが発生したとしても、外力のトルクさえ働かなければ角運動量は一定のままである、ということをこの式は示しています。このような法則を**"角運動量保存則"**といいます。
2.1.6 質量中心からみた運動量/角運動量.
質量中心の位置からみた、各質点の運動について考えてみましょう。座標系の原点から質量中心までのベクトルを＜Ｒ＞、質量中心から質点＜i＞までのベクトルを＜ｒi'＞として、質点の位置を＜ｒi＞は
［ 式2.20 ］ｒi＝Ｒ＋ｒi'
で表せます。さらに、質量中心の速度を＜Ｖ＞、質量中心からみた質点の相対速度を＜ｖi'＞として、質点の速度＜ｖi＞は
［ 式2.21 ］ｖi＝Ｖ＋ｖi'
で同様に表せます。式2.20の両辺に質点の質量＜ｍi＞を掛けてみると、
［ 式2.22 ］ｍiｒi＝ｍiＲ＋ｍiｒi'
となり、質点系の全質点の式2.22の総和をとると
［ 式2.23 ］Σｍiｒi＝ΣｍiＲ＋Σｍiｒi'
となります。ここで、＜ΣｍiＲ＞は合計質量＜ＭＲ＞と等しいので
［ 式2.24 ］Σｍiｒi＝ＭＲ＋Σｍiｒi'
さらに式2.9で左辺を置き換えると＜ＭＲ＞となり、結果として
［ 式2.25 ］Σｍiｒi'＝０
が導かれます。さらに、この式2.25の一階微分は
［ 式2.26 ］Σｍiｖi'＝０
となることも分かります。つまり、質量中心を基準とした場合の全質点の運動量の合計は必ず０になります。質量中心を原点とした座標系において、全質点の運動は平均すると"無いもの"となるわけです。これがどのような意味を持つかは、次節の剛体を学んだ時点で理解できることでしょう。次に質量中心を基準点にとった角運動量＜Ｌ＞を考えます。
［ 式2.27 ］Ｌ = Σ（ ｒi Ｘ ｍiｖi ）
= Σ{ （ ｒi' + Ｒ ） Ｘ ｍi（ ｖi' + Ｖ ） }
= Σ（ ｍiＲ ） Ｘ Ｖ + Σ（ ｍiｒi' ） Ｘ Ｖ + ＲＸΣ（ ｍiｖi' ） + Σ（ ｍiｒi' ） Ｘ ｖi'
この中で、＜Σ（ ｍiｒi' ）＞と＜Σ（ ｍiｖi' ）＞は式2.24と式2.25からそれぞれ０であることが分かっているので、消去して
［ 式2.28 ］Ｌ = Σ（ ｍiＲ ） Ｘ Ｖ + Σ（ ｍiｒi' ） Ｘ ｖi'
さらに＜Σ（ ｍiＲ ）＞は＜ＭＲ＞と等しいので、
［ 式2.29 ］Ｌ = ＭＲ Ｘ Ｖ + Σ（ ｍiｒi' ） Ｘ ｖi'
となります。このうち、右辺第一項は座標系の基準点から質量中心までの角運動量であり、右辺第二項は質量中心からみた各質点の角運動量の総計です。このように、角運動量は基準点から質量中心までの、いわば地球の公転を意味する値と、質量中心から各質点までの、いわば地球の自転を意味する値に分けられるのです。
2.2 剛体.
2.2.1 剛体とは.
さあ、ようやく剛体に話を進める準備が整いました。質点系などという不得要領な概念を持ちだして遠回りをしているかのように見えますが、そんなことは勿論ありません。なぜなら、剛体も質点系の特殊な一例に過ぎないからです。本章の始めに、剛体の定義を「大きさを持ち変形することがない固体」であるとの説明を行いました。これをもう少し正確な表現で記して、"任意の二質点間の距離が常に一定な質点系"である、とします。前節での質点系の話では、各質点毎の制約条件は特に定めていませんでした。対して剛体は、"相互距離を常に一定とする"という最も堅い制限を設けます。この制限によって、質点系に一意な形状、一意な姿勢を定めることができるようになるのです。例として、実際に頭の中で立方体を思い浮かべてみてください。立方体の８頂点にそれぞれ質点を一つずつ配置して、それぞれの相互距離を常に一定であると制限を設けました。このうち一つの質点をどんな風に動かしても、この制限がある限り立方体の形状が崩れることは絶対にありません。（ですよね？）このような"剛体"という概念を用いて、堅い物体の理想化を行います。
［ 図2.8 ］なんかの物体の二点を取りだして一定！な絵.
今までの話の流れから、剛体を「いくつかのボール同士を棒で繋いで接着剤で固定されている」、よく分子配列の模型にみられるようなイメージで考える方が多いかもしれません。これも間違いではありませんが、それよりも「一つの固体が原子や分子といったミクロなレベルで無数の粒子状に細分化されている。その粒子の一つ一つが質点であると考え、質点同士が強く結合している」といったイメージの方がよく使われます。どちらにしても定義に違いはありません。見え方に違いがあるだけです。
［ 図2.9 ］なんかの物体が細切れになってて、「質点として」「くっついているから剛体」な絵.
剛体の位置を正確に示すには、質点のような１点の位置だけでは不十分です。姿勢もきちんと把握するために、線形従属（※）でない三点の位置が最低限必要となります。（※） 線形従属とは、大雑把に言ってスカラー倍して同じ値になるベクトルを指す。座標（1,2,3）と（2,4,6）は線形従属であり、（2,2,3）は線形独立である。
- 一つの点のデカルト座標系での点Ａ（x,y,z : 情報数３）
- Ａの周りを極座標系で回る点Ｂ（s,t : 情報数２）
- ベクトルＡＢの周りを回る点Ｃ（q : 情報数１）
以上の計６情報数が必要条件であることから、"剛体の自由度は６である"などと呼ばれます。
［ 図2.10 ］剛体に一つの点と、それに極座標で回転してる点、円筒座標系で回転してる点、な絵.
2.2.2 並進運動と回転運動.
剛体は質点系の一種として定義された概念ですから、前節までに学んだ質点系の性質はそのまま適用されます。よって、剛体の質量中心の運動は、全質量が質量中心にある１質点と同様に扱う事が出来ます。さらに「二質点間の距離が一定である」という定義から、質量中心と質点との距離も常に一定です。しかも剛体内の全ての質点は、質量中心からみて全く同じ方向への回転運動を行います。これらのことから、剛体の運動は"剛体の質量中心の運動"と"質量中心からみた全質点に共通した回転運動"に分けて考えることができるのです。なお、前者は剛体の位置を決定づける運動であり、一般に"並進運動"と呼ばれます。対して後者は、剛体の姿勢を決定づける運動であり、一般に"回転運動"と呼ばれます。
［ 図2.11 ］剛体の運動は二つに分けろ！絵。
このように、剛体の運動を二つに分けて考える方法は、よくゲームなどで３Ｄ空間におけるローカル座標系からワールド座標系への変換行列を「平行移動（Translate）行列 x 回転（Rotate）行列 x 拡縮（Scale）行列」と分離して記述することとよく似ています。無論、剛体は大きさが変わらないと言う定義から、拡縮行列成分を必要としないのは自明ですね。
剛体の並進運動に関しては、一章で習ったとおりの質点の運動と見なすことで説明が付きます。ですから、これからの話から並進運動は除外します。残すは剛体の回転運動です。これも結局は並進運動と同様に、ニュートンの法則から導き出していくことになります。
2.2.3 角速度ベクトル.
回転運動に限って話を進める限り、質量中心は常に固定された不変な位置として扱うことができます。というのは、回転運動の計算が終わった後で並進運動の結果を全ての質点の位置に足してやれば、並進運動に関しての話は済んでしまうからです。よって、質量中心が固定された、回転しかしない物体の運動を考えればそれでよいのです。そこでまずは"回転"という情報を、どのようなパラメータで表現するべきか決める事にしましょう。
３Ｄプログラミングに慣れている方であれば、回転を表す方法として行列や四元数を用いることを真っ先に頭に浮かべるかと思います。しかし今回は、もっとシンプルにベクトルの形で表すことにします。まずは最初に、一つの定理を紹介しましょう。同形状をした、向きだけが異なる二つの剛体の姿勢を考えます。このうち片方の剛体に対して、質量中心を通る適切な回転軸を用意してやり、その軸周りの回転を一度行います。この一度の回転だけで、必ずもう片方と全く同じ方向を向かせることができるのです。これは「オイラーの定理」といい、どんな方向を向いた剛体同士であろうと成立します。この定理の証明は、x.x.x項にて後述します。
［ 図2.12 ］オイラーの定理な絵.
このオイラーの定理を用いて、ベクトルの形で回転を表現することが可能であることが分かります。この定理で必要としている情報は"回転軸の方向ベクトル"と"回転軸に対して回転する角度"ですから、回転角を大きさとする回転軸方向のベクトル＜φ＞を用意すればよいのです。剛体の初期状態である姿勢からどのような回転を行ったか、この回転角ベクトル＜φ＞を用いることで任意の姿勢を表現することが出来ます。なお、剛体の姿勢が変化することを考えると、この回転角ベクトル＜φ＞も時間によって変化する関数です。そこで姿勢の変化する速度を、一階微分の形で
［ 式2.30 ］ω = dφ/dt
と書くとします。この＜ω＞は、剛体の回転角ベクトル＜φ＞の時間変化率を表しています。並進運動での位置ベクトルの時間変化率を速度ベクトルと呼んだように、＜ω＞は"角速度ベクトル"と呼ぶこととします。同様に回転角ベクトルの二階微分
［ 式2.31 ］β = dω/dt = d^2φ/dt^2
の＜β＞を"角加速度ベクトル"と呼ぶのが自然でしょう。回転角、角速度、角加速度の関係は、並進運動での位置、速度、加速度の関係と相似であると見なせます。角速度ベクトル、角加速度ベクトルが求まれば、現在の剛体の回転の様子が求まります。つまるところ剛体の回転運動を追いかけることとは、これらの値を求めることに他なりません。
ここで、角速度ベクトル＜ω＞と剛体内の各質点との関係について触れておきます。
［ 図2.13 ］角速度ベクトルと質点の速度の対応絵.
剛体内の適当な質点をとり、質量中心からの位置ベクトルを＜ｒ＞とします。この位置ベクトル＜ｒ＞は、回転軸＜ω＞の周りを半径＜|ｒ|sinφ＞で円運動の軌跡を描きます。さらに、並進運動が静止状態であると仮定した場合の速度ベクトル＜ｖ＞は、角速度ベクトル＜ω＞に対して垂直で、位置ベクトル＜ｒ＞の位置から角速度ベクトル＜ω＞へ引いた方向ベクトルに対しても垂直であることが図2.13から見て取れます。また、速度ベクトル＜ｖ＞の大きさ＜|ｖ|＞は＜|ｒ||ω|sinφ＞です。これらにより、質点の位置ベクトル＜ｒ＞と角速度ベクトル＜ω＞の関係は外積の形で記述することができ、最終的に以下の式が成り立ちます。
［ 式2.32 ］dｒ/dt = ｖ = ωＸｒ
この式を用いることで、角速度ベクトルから任意の質点の速度ベクトルを速やかに決定的に求めることが出来ます。これは後々まで利用される重要な式となります。
2.2.4 角速度ベクトルの反映.
角速度ベクトルから回転角ベクトルの変化を求める方法は、速度ベクトルから位置ベクトルの変化を求めるのと全く同様の手順、やはり微分方程式を解くことになります。微分方程式の解法は一章で習ったとおりです。特に新しい手法が必要なわけではありません。しかし実際に３Ｄプログラミングを行う場合において、回転を表現する手法として回転角ベクトルはあまり適しているとは言えません。それよりも行列表現や四元数の方が便利ですし、何より普段から使い慣れていることでしょう。ですから、角速度ベクトルからは回転角ベクトルの変化を求めるのではなく、回転用の行列や四元数の変化を求める方法を探っていくべきです。
まずは回転行列の変化について説明します。3行3列の回転行列＜Ｒ（t）＞は、本質的には剛体の持つ座標系の基底ベクトル＜ｅi（t）＞を成分毎に３つ並べただけのものと考えられます。この基底ベクトルは単位化されており、それぞれ直交しているとします。
［ 式2.33 ］ｅx（t） = ［ｒ00 ｒ01 ｒ02］
ｅy（t） = ［ｒ10 ｒ11 ｒ12］
ｅz（t） = ［ｒ20 ｒ21 ｒ22］
［ 式2.34 ］Ｒ（t） = ［ｒ00 ｒ01 ｒ02］ = ［ ｅx（t） ］
［ｒ10 ｒ11 ｒ12］ ［ ｅy（t） ］
［ｒ20 ｒ21 ｒ22］ ［ ｅz（t） ］
ここで、基底ベクトル＜ｅi（t）＞の一階微分は式2.32を使うことで求めることができ、回転行列＜Ｒ（t）＞の一階微分は
［ 式2.35 ］dＲ（t）/dt = ［ dｅx（t）/dt ］ = ［ ω（t）Ｘｅx（t） ］
［ dｅy（t）/dt ］ = ［ ω（t）Ｘｅy（t） ］
［ dｅz（t）/dt ］ = ［ ω（t）Ｘｅz（t） ］
となります。このように基底ベクトルの集合として考えてやることで、行列の変化をシンプルに表すことができるのです。なお、今回のようなベクトルと行列の外積は、行列と行列の積に変形させることで、見た目に綺麗にすることが可能です。が、代わりに計算量は増えてしまいますので、よほど行列の積が高速なハードウェア上で動かすのでない限り、忘れてしまって全く構いません。
このような、角速度ベクトルを回転行列に反映させる方式には弱点があります。数値解法で逐次計算していくことで、各基底ベクトルが計算を続けていく度に累積誤差を抱えてしまうのです。これにより、基底ベクトルが単位ベクトルではなくなってしまったり、直交しなくなってしまったりしてしまいます。誤差を隠すためには、行列の正規化を行ってやる必要があります。正規化は毎回行う必要はありませんが、どの程度の誤差を許してよいか判断が難しいところです。
回転行列ではなく回転四元数に反映させる手法では、この誤差を考えずに済みます。角速度ベクトル＜ω（t）＞から回転四元数の変化分に変換する式は以下の通りです。
［ 式2.36 ］ｑ（Δt） = ＜1,1/2Δtω（t）＞
さらに、現在の回転四元数からの更新は
［ 式2.37 ］ｑ（t+Δt） = ｑ（t）ｑ（Δt）
という式で表すことが出来ます。回転四元数のままでは平行移動行列などとの結合ができないので、ここからさらに回転行列の形に戻してやる必要があります。とはいえ、計算誤差の問題を考慮に入れる必要が無いという利点がありますので、何かしらの懸念事項がなければ極力四元数を利用するべきでしょう。
2.2.5 慣性テンソル.
前項にて角速度から剛体の姿勢を求める方法を説明しましたので、次に角運動量と角速度との関係について考えていくことにしましょう。ここからの話、角運動量と書いた場合には全て"質量中心を基準点とした角運動量"のことを暗黙的に指します。並進運動での運動の勢いを運動量と呼ぶように、回転運動での運動の勢いを表した値が角運動量なのです。質点系の全角運動量＜Ｌ＞は、式2.15のまま
［ 式2.38 ］Ｌ = Σ（ ｒi Ｘ ｐi ）
= Σ（ ｒi Ｘ ｍiｖi ）
でした。剛体も質点系の一種ですから、この式は論無くそのまま利用できます。式2.38の＜ｖi＞は質点＜i＞の速度ベクトルであり、これは角速度ベクトル＜ω＞を用いて
［ 式2.39 ］ｖi ＝ ω Ｘ ｒi
と書けることは2.2.3項にて導きました。式2.38と式2.39を組み合わせると、
［ 式2.40 ］Ｌ = Σ（ ｒi Ｘ ｍi（ω Ｘ ｒi） ）
となります。ここで、ベクトルの三重積の公式
［ 式2.41 ］（ａ Ｘ ｂ）Ｘ ｃ ＝ （ａ・ｃ）ｂ－（ｂ・ｃ）ａ
を用いて、
［ 式2.42 ］Ｌ = ω（Σ（ ｍi|ｒi|^2））－Σ（ ｍiｒi（ｒi・ω））
と変形させます。この＜Ｌ＞はベクトルなので、これをxyzの成分表示にしてみると
［ 式2.43 ］Ｌx ＝ ＋ ωxΣｍi（ｙi^2 + ｚi^2） － ωyΣｍiｘiｙi － ωzΣｍiｘiｚi
Ｌy ＝ － ωxΣｍiｙiｚi ＋ ωyΣｍi（ｘi^2 + ｘi^2） － ωxΣｍiｙiｘi
Ｌz ＝ － ωzΣｍiｚiｘi － ωyΣｍiｚiｙi ＋ ωzΣｍi（ｚi^2 + ｙi^2）
となります。なんとも見通しが悪いので、それぞれの角速度に掛かっている係数を以下のような記号＜Ｉ**＞に置き換えます。
［ 式2.44 ］Ｉxx ＝ ＋ Σｍi（ｙi^2 + ｚi^2）
Ｉxy ＝ － Σｍiｘiｙi
Ｉxz ＝ － Σｍiｘiｚi
Ｉyx ＝ － Σｍiｙiｘi
Ｉyy ＝ ＋ Σｍi（ｘi^2 + ｘi^2）
Ｉyz ＝ － Σｍiｙiｚi
Ｉzx ＝ － Σｍiｚiｘi
Ｉzy ＝ － Σｍiｚiｙi
Ｉzz ＝ ＋ Σｍi（ｚi^2 + ｙi^2）
これらを用いて式2.43を書き直すと、
［ 式2.45 ］Ｌx ＝ ωxＩxx + ωyＩxy + ωzＩxz
Ｌy ＝ ωxＩyx + ωyＩyy + ωzＩyz
Ｌz ＝ ωxＩzx + ωyＩzy + ωzＩzz
と幾分すっきりと書けました。さらに、次のような行列＜Ｉ＞を考えます。
［ 式2.46 ］Ｉ = ［ Ｉxx Ｉxy Ｉxz ］
［ Ｉyx Ｉyy Ｉyz ］
［ Ｉzx Ｉzy Ｉzz ］
この行列を用いると、式2.45は
［ 式2.47 ］Ｌ = Ｉ・ω
と、非常にシンプルな形で記述することができました。この行列＜Ｉ＞を"慣性テンソル"といいます。加えて慣性テンソルの対角成分（Ixx,Iyy,Izz）を慣性モーメント、その他の成分（Ixy,Ixz,Iyx,Iyz,Izx,Izy）を慣性乗積といいます。この式2.47が示すとおり、慣性テンソルは"角速度ベクトルと角運動量ベクトルの関係を結ぶパラメータ"です。この関係は、並進運動でいうところの運動量を求める式1.13、
［ 式1.13 ］ｐ = ｍｖ
に非常によく似ています。並進運動での質量＜ｍ＞が"位置の変化しづらさ"を表しているように、回転運動での慣性テンソル＜Ｉ＞は"回転の変化しづらさ"を表しています。トルクは角運動量の時間変化率であり、これは運動量の時間変化率が力であることとアナロジーの関係にあります。角運動量から角速度ベクトルを求めるのに慣性テンソルが関係を保つのも、運動量から速度ベクトルを求めるのに質量が関係を保つのもアナロジー的です。並進運動では力から位置の変化を求めますが、回転運動ではトルクから姿勢の変化を求めます。結局のところ回転運動は、式は少し複雑になったもののやっていることは並進運動と大して変わっていない、という点は見落とせないポイントでしょう。角運動量と角速度ベクトルの関係性が分かれば、あとは前章と同様に、微分方程式を解析解法しかり数値解法しかりの方法で解いて、トルクから角速度ベクトルを求めてやればよいのです。
回転角ベクトル φ ＜-＞ 位置ベクトル ｒ
角速度ベクトル ω ＜-＞ 速度ベクトル ｖ
角加速度ベクトル β ＜-＞ 加速度ベクトル ａ
慣性テンソル Ｉ ＜-＞ 質量 ｍ
角運動量 Ｌ ＜-＞ 運動量 ｐ
トルク Τ ＜-＞ 力 ｆ
式 Ｉω = Ｌ ＜-＞ 式 ｍｖ = ｐ
式 dＬ/dt = Τ ＜-＞ 式 dｐ/dt = ｆ
しかし、剛体に関する話はこんなに簡単に片づくわけではありません。実はまだ厄介な問題が一つ残っています。それは慣性テンソル＜Ｉ＞の中身です。並進運動での運動の変化しにくさを表す質量＜ｍ＞は、不変なスカラー値として扱っても問題はありませんでした。しかし、困ったことに慣性テンソル＜Ｉ＞は式2.44を見れば分かるとおり、剛体内の質点の位置がテンソルの中に埋め込まれてしまっています。すなわち、時間が経って姿勢が変化したとしたら、当然質点の位置も変化します。結果として、慣性テンソルの中身も変化してしまうのです。毎回姿勢が変わるたびに慣性テンソルを計算していては、計算量の大きな無駄になります。というわけで、この問題の解決法を探っていきましょう。
［ 図2.14 ］慣性テンソルが回転すると変わってしまうよ！絵.
2.2.6 慣性主軸と対角化.
前述した慣性テンソルをよくよくみると、慣性乗積の部分が全て Ijk = Ikj となっているのが分かります。このようなシンメトリーな行列を"対称行列"といいます（全ての値が実数なので"実対称行列"ともいいます）。そして対称行列は"対角化"が可能です。対角化とは、行列Ｍに対してＡ^-1ＭＡが"対角行列"となるようなＡを見つけることをいいます。対角行列とは、j=kでないＭjkが全て0な行列を指します。慣性テンソルで言えば、慣性乗積の部分が全部0であればそれは対角行列と呼べるわけです。
あまり耳慣れない言葉がずらりと並びましたが、つまりは「慣性テンソルをうまく加工すると、慣性乗積が全て0になる新しい行列を作ることが出来る」と言いたいのです。このうち、うまく加工すると、というのを詳しく説明すると、「慣性テンソルにとって都合の良い座標系をもう一つ作る」ことなのです。これは３Ｄモデルの頂点位置を指示するのに、ワールド座標系ではなくローカル座標系を使用するのと同じです。この、慣性テンソルが対角化できる座標系のことを"慣性主軸座標系"と呼びます。慣性主軸座標系内では慣性乗積が0なのです。2.8項の長ったらしい式もずいぶんスッキリしそうですよね。この慣性主軸座標系内での慣性モーメントを"主慣性モーメント"と呼びます。慣性乗積分の情報が慣性モーメントに集まってより重要な情報になったわけです。
さぁ、ここまで来て実際の慣性テンソルを対角化方法は...というと、大学の線形代数で勉強する、実対称行列の固有値問題を解くことになります。解析的に解けるものは殆ど無く、大抵は反復的に解く必要があり、深く解説してしまうと話がややこしくなりすぎてしまいますのでここでは割愛します。
2.2.7 オイラーの運動方程式.
慣性主軸座標系での剛体の回転運動を追いかけていきます。まず、どんなものでも構わないので適当なベクトル＜Ａ＞があるとします。このベクトル＜Ａ＞が剛体と一緒に回転しているとすると、通常の（ワールド）座標系から見た一階微分は式2.32から
［ 式2.48 ］dＡ/dt = ωＸＡ
となります。さらに、ベクトル＜Ａ＞が剛体からみて変化しているとします。この場合、式2.48に時間変化率＜δＡ/δt＞を足してやることになるので、
［ 式2.49 ］dＡ/dt = ωＸＡ + δＡ/δt
となります。これはベクトル＜Ａ＞の運動を、剛体に対する追従（右辺第一項）と剛体に対する変化（右辺第二項）に分離して考えた形になっています。この式2.49を用いて、剛体の角運動量＜Ｌ＞の時間変化率＜dＬ/dt＞、つまりトルク＜Τ＞を表すと以下のようになります。
［ 式2.50 ］dＬ/dt = Τ = ωＸＬ + δＬ/δt
さらに式2.47を用いて
［ 式2.51 ］Τ = ωＸＬ + δＩ・ω/δt
と変形できます。慣性主軸座標系では慣性テンソルは不変であるとしてよいのですから、
［ 式2.52 ］Τ = ωＸＬ + Ｉδω/δt
という式が導かれます。＜δω/δt＞を角加速度＜β＞で置き換え、xyzの成分表示にしてみると以下のような式になります。このとき慣性テンソル＜Ｉ＞は慣性主軸座標系であることから、慣性乗積は０としています。
［ 式2.53 ］Τx = （ Izz - Iyy ） ωzωy + Ixx・βx
Τy = （ Ixx - Izz ） ωxωz + Iyy・βy
Τz = （ Iyy - Ixx ） ωyωx + Izz・βz
この式2.53を"オイラーの運動方程式"といいます。求めたいのは角速度ベクトルの時間変化率＜β = δω/δt＞なので、もう少し式を変形して
［ 式2.54 ］βx = （ Τx + （ Iyy - Izz ） ωzωy ） / Ixx
βy = （ Τy + （ Izz - Ixx ） ωxωz ） / Iyy
βz = （ Τz + （ Ixx - Iyy ） ωyωx ） / Izz
となります。これでトルク＜Τ＞から直接、角加速度＜β＞を求めることが出来るようになりました。ただしこの角加速度＜β＞は慣性主軸座標系での値ですから、後で通常の（ワールド）座標系に変換してやる必要はあります。角加速度から角速度の求め方はいつものように数値解法に頼ることになるでしょう。
このようにオイラーの運動方程式を用いることによって、慣性テンソルの再計算と言った煩雑な計算をせずに回転運動を計算することができるようになります。これで与えられた力からトルクが生まれ、トルクから角速度を求め、現在の姿勢からどのような変化が起きるか、つまりは"回転運動"を完全に把握することが出来るようになったのです。
2.2.8 独楽の歳差運動.
本章の最後に、独楽の歳差運動について触れておきましょう。回転していない素な状態の独楽は、立たせたまま置こうとしてもすぐに倒れてしまいます。にも関わらず、なぜか回転している独楽はふらついたような不安定な動きをしつつも直立の体勢を崩しません。普通に考えたら倒れてしまいそうな状態でも、回転していれば倒れずに踏ん張るのです。歳差運動とは、このような回転している独楽の軸部分に見られる、円錐を模した動きをいいます。この不可思議な運動はどのようにして起こるのでしょうか？
［ 図2.15 ］独楽が回転して歳差運動を起こしている絵.
これは、角速度ベクトルの合成で説明が付きます。現在回転している独楽の勢いを角速度ベクトルで表すと、当然の事ながら独楽の軸に沿ったベクトルとして表現されます。直立していれば鉛直方向、傾いていればその傾き方向に対して角速度ベクトルは向きを変えます。このときの重力の働きを考えます。直立した状態であればそのまま真下方向に対して力が働き、独楽の地面との接触点ではトルクを得ることはありません。よって直立した状態であれば回転に対して倒れる方向への回転の勢いは生まれません。では、少し独楽が傾いたときにはどうなるでしょうか？重力によって、傾いた方向へのトルクが働きます。このトルクの発生する位置は独楽の地面との接触点であり、このトルクのベクトルは傾いている方向の"逆方向"に対して働きます。このトルクによって、もし回転していなければ傾いている方向への回転が行われ最終的に倒れてしまいます。しかし独楽が回転しているのであれば、前述した軸に沿った角速度ベクトルとの合成が行われ、角速度ベクトルは最終的に鉛直方向へ向きを変えます。この角速度ベクトルを得た独楽の軸は、重力とは垂直の方向に動きを進め、鉛直方向への軸を中心としたゆっくりとした回転をするわけです。
［ 図2.16 ］図2.15にベクトルを足してベクトルの合成で鉛直方向に太いベクトルを描いた絵.
以上の話から、いくら独楽であっても空中に浮いた状態で歳差運動を行うことはできないことが分かります。重力だけではなく、地面からの反発力がなければ歳差運動は行われないのです。これで独楽の歳差運動についての説明は付きました。が、それではこの"重力"とか"反発力"などといった力はどこからどのように生まれるのでしょうか？
次章は「力」について探っていきましょう。
終わりに
長すぎるので序盤戦も序盤戦のみ、となってしまいました。この続きはまたどこかで...
それではまた来年、12/24のアドカレでお会いしましょう...
［eof］