3.壁との衝突

壁との弾性衝突によるボールの速度変化は、衝突公式(16)で与えられる。
この章のシミュレーション(クリックで計算開始)

弾性衝突による速度の変化を知りたい

壁面との衝突を含む場合のボールの運動 x(t) が知りたい(右図)。ボールは、衝突前、重力などの力を受けて自由に運動する。その後、壁面に衝突して跳ね返る。跳ね返った後は、次の衝突まで再び自由に運動する。衝突の瞬間以外は、第2章の場合そのものなので、すでに計算方法は分かっている。従って、新たに考える必要があるのは、衝突時の運動の変化のみである。特に、衝突直後の速度 ˙xout が分かれば十分である。なぜなら、運動を一意的に決めるのに必要な初期値はボールの位置・速度であるため(第1章の1.1節)、衝突直後の位置・速度が分かれば、それらを初期値としてその後の運動が一意的に決まるからである(未知なのは速度のみ)。
ただし、ボールの弾み方は、ボールや壁の素材によって異なり、よく跳ねるものもあれば、そうでないものもある。そのため、「弾みやすさ」を決めてやらないと ˙xout は決まらないはずである。ここでは、最もよく弾む理想的な場合を考えることにする。これを弾性衝突という。(このような理想的な極限が存在することは直感的にわかるだろう。例えば、床にボールをそっと落とした時、元の高さまで戻ってくるのが弾性衝突である。)
弾性衝突での速度の変化は、直感的に分かる。即ち、右図のように、衝突前の速度 ˙xin に対し、「壁と垂直な方向の成分を反転」したものが ˙xout になる(実験的にも確認できる)。従って、これを認めてしまえば、後は「…」部分をどのように求めるかという、幾何学的な問題になる。そのようにして、衝突公式を得ることができるはずだが、それを実際に使うには、衝突点での壁面の向きなどを与える必要がある。そのためには、壁をどのように数式で表すかを考える必要がある。
そこでこの章では、 x(t) の計算を、以下の3つの節に分けて議論を行う: 3.1:3.2:3.3:x(t)

3.1 衝突公式(2)の導出

この節では、まず幾何学的な考察により、衝突後の速度 ˙xout を与える衝突公式(2)を導く。その後、力学的な考察から得られるダランベールの原理(5)と弾性衝突条件(6)を用いて、同じ衝突公式を再導出する。

衝突後の速度 ˙xout は、「壁と垂直な方向の成分」が反転する:衝突公式(2)

平らな硬い床に、よく弾むボールを垂直に落としてみる。すると、経験的に分かるように、衝突の前後でボールの運動方向は反転するが、速さは変わらない。一方、斜めにぶつけた場合、衝突後の速度 ˙xout のうち、「壁と平行な成分」は衝突前と変わらず、「床と垂直な成分」は反転する。(現実的には、平行成分も摩擦の影響で多少変化するが、ここでは無視できるとする。)
これを数式で表すためには、衝突前の速度 ˙xin を、床や壁に対して「平行な成分 ˙xin 」と「垂直な成分 ˙xin 」に分解すればよい: ˙xin=˙xin+˙xin そうすれば、衝突により垂直成分 ˙xin が反転するのだから、 ˙xout は以下のようになる: ˙xout=˙xin˙xin(1)
各成分 ˙xin,˙xin を求めたいわけだが、これは、「壁面への正射影行列 P 」および「壁面に垂直な方向への正射影行列 P 」を用いて、それぞれ ˙xin=P˙xin,˙xin=P˙xin と書ける(以下の【3.1-注1】の式(3))。これらを上式(1)に代入すれば ˙xout=(PP)˙xin=(12ˆnˆnT)˙xin(4)(2) となる。ただし、 ˆn は、壁に垂直な単位ベクトル(=単位法線ベクトル)である。よって後は、衝突点での ˆn を求めることさえできれば、上式(2)が確定し、 ˙xout が計算できることになる。

【3.1-注1】正射影行列:式(4)

右図のように、任意のベクトル v を、ある平面 S に対して「平行な成分 v 」と「垂直な成分 v 」の和に分解する: v=v+v この時、各成分 v,v は、 S と垂直な単位ベクトル ˆn を用いて、それぞれ以下のように書ける: v=(1ˆnˆnT)Pv,v=(ˆnˆnT)Pv(3)
1 は単位行列、 ˆnˆnT は以下の【3.1-注2】の外積である。)上式に現れた行列 P,P のことを、正射影行列という[1]P=1P,P=ˆnˆnT(4)

補足

[1] 任意のベクトル v に対し、 Pv は「 v を平面 S に正射影」したもの、 Pv は「 S に垂直な方向に正射影」したものになるわけである。例えば、 Sx-y 平面である場合、 nz 軸方向の単位ベクトルを取ればよいので、正射影行列(4)は以下のようになる: P=100010000,P=000000001 Px-y 平面への正射影行列(= x,y 成分のみを取り出す行列)Pz 軸への正射影行列(= z 成分のみを取り出す行列)となっている。

【3.1-注2】内積と外積

2つのベクトル a,b に対し、 abT は行列の積として定義できる:bTb の転置) abT=axayaz[bxbybz]=axbxaxbyaxbzaybxaybyaybzazbxazbyazbz 転置記号 T が外側にあるので、これを外積(outer product)という。
逆に、 T が内側にある aTb を考えることもでき、内積(inner product)という: aTb=[axayaz]bxbybz=axbx+ayby+azbz

補足

内積は ab と表されることもあるが、本サイトでは使わないことにする(「 」記号は積を区切るためだけに使う)
外積の例として、式(4)の正射影行列 P を考えてみる。 P に任意のベクトル v を掛けると、以下のようになる: Pv=(ˆnˆnT)v=ˆn(ˆnTv)=(ˆnTv)ˆn (1行目の変形は行列の積の結合律を用いた。)最後の式は、大きさが ˆnTv で、向きが ˆn であることを示している。 ˆnTv は、内積の性質により、 vˆn 方向成分なので、確かに ˆn 方向の正射影になっていることが分かる。

力学的条件:ダランベールの原理(5)と弾性衝突条件(6)

衝突公式(16)はすでに得られているが、幾何学的すぎるので、力学的に再導出することを考える。前章でみたように、速度の変化を知るためには加速度が分かればよく、加速度を知るためにはボールに働く力が分かればよい。とはいえ今の場合、速度が瞬間的に変化するので、速度の時間微分である加速度は定義できないし(第1章の【1.1-注1】補足参照)、衝突時の力を完全に決定することも容易ではない。そのため、速度の変化を求めるにあたって、これまでの考え方が使えない。
壁に衝突したボールが速度 ˙x の向きを変えるのは、壁にめり込むのを防ぐために働く拘束力 fc のためである(添え字 c は拘束を意味するconstraintの頭文字)。上述の通り、衝突の前後で変化するのは ˙x の「壁に垂直な成分」のみである。速度が変化する方向と力の向きは一致するので、 fc の向きが壁に垂直となる。即ち、 fc は以下のように書ける:λ は未定乗数と呼ばれる未知数) fc=λn(5) これをダランベールの原理と呼ぶ。未定乗数 λ は、衝突の瞬間にのみ 0 でない値を持つが、衝突の瞬間の短い間にも急激に変化するため、直接測定することは難しい。
ダランベールの原理(5)は、弾性衝突でなくても成り立つので、これとは別に、弾性衝突であるための条件があるはずである。それは、衝突の前後でボールの運動の勢いが減衰しないこと、即ち、速度の大きさが不変: |˙xout|=|˙xin|(6) であるという条件である。これを弾性衝突条件という。

衝突公式(16)の力学的な再導出

示したいのは、ダランベールの原理(5)と弾性衝突条件(6)から、弾性衝突の公式(16)が導かれることである。そのために、ダランベールの原理(5)を、衝突前後の速度 ˙xin,˙xout に対する関係式にしたい。ダランベールの原理により、拘束力の向きは壁に垂直なのだから、速度の変化も壁に垂直である。よって、以下の関係式が成り立つ:Λ は未知数) ˙xout=˙xin+Λn(7)
後は、未知数 Λ を求めればよいわけだが、そのためには弾性衝突条件(6)(を2乗したもの)˙xTout˙xout=˙xTin˙xin(8) を用いればよい。具体的には、式(7)を式(8)に代入して ˙xout を消去すればよい。すると Λ の単純な2次方程式になり、その解は以下の2つである: Λ={02|n|2nT˙xin Λ=0 の場合、速度が変化しないため、壁にめり込んで拘束条件を破ってしまう。よってこれは不適当なので、もう1つの解が正しい解である。この Λ を式(7)に代入すると、 ˙xout が確定し、幾何学的に求めた式(16)に一致することが分かる。

3.2 壁の定式化

この節では、壁の形状の表現方法について述べた後、式(2)で計算方法がまだ分かっていなかった法線ベクトル n の計算方法を与える。

拘束条件(11)によって壁の形状を表す

壁面の法線ベクトル n を求めたいのだが、もちろんこれを計算するには、壁の形状が予め与えられている必要がある。そのために、壁の形状を、関数 G(x) で表すことにする。 G(x) は、ボールの位置 x が壁の外にある時に正の値を取り、壁の内部で負の値を取る(右図)。つまり、 x は、不等式 G(x)0(9) を満たす点のみを取ることができ、不等式が破れる点が壁になる。この不等式、または G を拘束条件と呼ぶ。例えば、半径 R の球の内部にボールを閉じ込める場合、ボール位置 x の取りうる範囲は |x|R で表されるので、拘束条件は G(x)R|x|0(10) とすればよい。ただし、 G の取り方には任意性があり、例えば GR2|x|2 などとしてもよい。
記号 は以前にも登場しているが、定義または恒等式であることを表すものである。特に、方程式でないことを強調する。例えば式(10)の左側は拘束条件 G の定義であり、 x に対する方程式ではない(右側の不等式は x に対する方程式である)。ただし自明な場合には、定義や恒等式であっても = と書く。

壁面に垂直なベクトル n :式(11)

後は、拘束条件 G(x) から、壁面の法線ベクトル n を求めればよい。これは純粋に幾何学的な問題であり、その答えは、以下の【3.2-注3】の式(12)で与えられる: n=GxGyGzG(11)
なお、 G0 になる可能性がある。例えば、任意の拘束条件 G に対し、 GG3 を考えてみる。 G は、壁の内部で負、外部で正となっているので、拘束条件として正しいが、壁面上で G=3G2G=0 となる( G=0 )。こうなるのは、 G の値が壁面上で停留値をとるような場合である。こうなってしまうと、法線ベクトルではなくなってしまう。上述のように G の取り方には任意性があるので、それを利用して G0 が成り立つようにしておく必要がある。

【3.2-注3】曲面に垂直なベクトル

ある曲面が G(x)=0 によって表されているとする。この時、曲面上の点 x における以下のベクトル GGxGyGzG(12) は、 0 でなければ、曲面に垂直である。 はナブラと読む。 xGGx の略記である。

導出

曲面上の2点 x,x+δx を考える(両方とも G=0 を満たす)δx を十分小さく取っておき、 G(x+δx) を1次近似すると、以下の【3.2-注4】の式(13)により、以下のようになる:(右辺の量は x での値) G(x+δx)0G0+dGdxδxdGdxδx0 接平面上の任意の δx に対してこの式が成り立つので、赤字部分の転置行列をとったもの: G(dGdx)T は、接平面と垂直であるか、あるいは 0 である。

【3.2-注4】多変数関数の1次近似:式(13)

右図のように、2変数 x=(x,y) の任意の関数 G=G(x) のグラフを考える。このグラフ上の任意の1点 x の近傍を拡大していくと、右図のように次第に平面(=接平面)に近づいていく。(これは、地球の表面は球面であるにもかかわらず、地上の我々からは平面に見えることからも分かる。)
平面のグラフは一般に1次関数で表せるので、 x の近傍の値 G(x+δx) も、微小なベクトル δx の1次関数で近似できる。実際、以下のようになる[1](右辺の量は x での値) G(x+δx)G+Gxδx+Gyδy=G+dGdxδxdGdx[GxGy](13)(14)
これは、1変数関数のグラフを1次近似したものが、接線になるのと同様である(第1章の【1.1-注1】)。
より多変数の場合にも自然に拡張できる。例えば、3変数の場合、式(14)を以下で置き換えればよい: dGdx[GxGyGz]

補足

[1] 式(13)の緑字部分がこうなることは、 δx=0 の時に両辺が等しくなることから分かる。赤字部分が式(14)のようになることについては、例えば、式(13)において、 δx=[δx0] としてみればよい。これは、 y を固定して、 Gx の1変数関数とみなすことに対応する。よって、 δx の係数は、 Gx (= x による偏微分)でなければならない。 δy の係数についても同様である。
式(13)のように1次近似ができる時、 G は、その点 x において微分可能(あるいは全微分可能)であるという。1変数関数の場合と同様、角ばった点があったりすると、その点では微分できない。
式(13)は、 δGG(x+δx)G(x) と定義することにすれば、シンプルに δGdGdxδx(15) とも書ける。微分の計算をする際にはこちらのほうが便利である。例えば、次節の【3.3-注1】。
なお、気にする必要はまずないが、 dGdx の各成分が存在しても(=偏微分が存在しても)、 G が微分可能であるとは限らない。というのも、偏微分が存在しているだけでは x,y 方向に1次近似できることしか言えず、別の方向では1次近似できないかもしれないからである。

衝突後の速度 ˙xout :式(16)

以上により、衝突後の速度 ˙xout の公式を得るには、衝突公式(2)に、法線ベクトル n の式(11)を代入すればよい。すると、以下のようになる:ˆGG を正規化したもの: ˆG=G/|G| ˙xout=[12ˆG(ˆG)T]˙xin(16) これが求めたかった公式である。これでようやく、 ˙xout が計算できるようになった。

3.3 ボールの運動 x(t) の計算

この節では、まず、計算方法についてまとめ、その後、いくつかの例に対して具体的な計算を行う。

衝突を含む運動 x(t) の計算方法

さて、もともと知りたかったのは、衝突を含むボールの運動 x(t) である。実際に x(t) を計算するには、まず、衝突する瞬間、即ち G(x(t))=0 となる瞬間の時刻 t まで、前章で述べた運動方程式 m¨x=f を用いて計算を進めるf は例えば重力)。壁に接している状態で、式(16):(再掲) ˙xout=[12ˆG(ˆG)T]˙xin(17) 上式を使って衝突後の速度 ˙xout を計算する。それ以降は、衝突位置と ˙xout を初期値として、再び衝突するまで m¨x=f を用いて計算を進めていけばよい。
なお、初期値の取り方によっては、壁に接触したまま、壁に沿うように運動する場合も考えられる。その場合は、衝突が一瞬で終わるという暗黙の仮定が成立しないので、これまでの議論が使えない。実際、式(17)をそのまま使うと、 ˙xin の壁面と垂直な方向成分は 0 なので、 ˙xout=˙xin となってしまい、壁が無いのと同じになってしまう。詳しくは後の第6章で扱うが、このような場合、拘束力を運動方程式に加えるという形で対応できる。

球・円と衝突する場合:式(18)

半径 R の円または球の内部にボールを閉じ込めた場合を考える。拘束条件 G およびその微分 G は、以下のようになる:(以下の【3.3-注1】の公式(19)を使う) GR|x|0G=dGdxT=ˆx これを式(17)に代入すれば、衝突後の速度 ˙xout を求める公式が得られる: ˙xout=(12ˆxˆxT)˙xin(18)
これを用いて数値計算を行うと、右図のようになる。円・球の半径 R はともに R=1m である。

【3.3-注1】 |x| の微分公式

位置ベクトル x の絶対値 |x| の微分は、以下のようになる: ddx|x|=ˆxT(19)

導出

微分の定義(15)を使えばよい。即ち、求める微分は、 δ|x| の1次近似: δ|x|=δyy|x|2dydyδy11δfdfdyδy=12|x|δydydy=12y=12|x|ˆxTδx⎪ ⎪ ⎪ ⎪⎪ ⎪ ⎪ ⎪δy=|x+δx|2|x|2=(x+δx)T(x+δx)xTx=2xTδx+δxTδx2xTδx+0 における δx係数である。

y=f(x) の形の床との衝突の場合:式(20)

もう1つの例として、2次元平面上において、 y=f(x) のグラフの形の床との衝突を考える(右図)。拘束条件 G およびその微分 G は、以下のようになる:fx での微分: f=df/dx G(x,y)yf(x)0G=[xGyG]=[f1] これを式(17)に代入すれば ˙xout=(12f2+1[f2ff1])˙xin(20) となる。
例えば、 f(x) として以下のもの:a,b,c は定数) f(x)=acosbx+cx2 をとることにして、数値計算を行うと右図のようになる。グリッドの間隔は 1m である。
ところで、衝突は計算誤差を増大させる性質を持つ。というのも、衝突により、位置のずれが運動方向のずれを生じさせるからである(床が湾曲している場合)。そのため、すぐに誤差が無視できなくなり、計算が意味をなさなくなってしまう。実際、上のシミュレーションでも、(ステップ幅 δt が毎回異なる実装のため)時間が経つにつれて毎回全く異なる運動になる。結果が異なってしまうのは理論や計算の欠陥のようにも見えるが、そうではなく、現実世界での再現困難性を反映しているのである。実際、曲がった床にボールを落として跳ねさせる実験を行うと、同じところから落としたつもりでも、数回跳ねると毎回違う運動になってしまう。この性質が、理論にも正しく反映されているわけである。