\newcommand{\sbr}[1]{\left[#1\right]} \newcommand{\br}[1]{\left(#1\right)} \newcommand{\cbr}[1]{\left\{#1\right\}} \newcommand{\abr}[1]{\left|#1\right|} \newcommand{\exbr}[1]{\left\langle #1 \right\rangle} \newcommand{\nbr}[1]{\left\lVert #1 \right\rVert} \newcommand{\sC}{\mathbb{C}} \newcommand{\sN}{\mathbb{N}} \newcommand{\sR}{\mathbb{R}} \newcommand{\sL}{\mathit{\Lambda}} \newcommand{\sSS}{\mathit{\Omega}} \newcommand{\sFW}{\mathfrak{W}} \newcommand{\sFB}{\mathfrak{B}} \newcommand{\sFF}{\mathfrak{F}} \newcommand{\sFP}{\mathfrak{P}} \newcommand{\Open}[1]{\mathfrak{O}(#1)} \newcommand{\Close}[1]{\mathfrak{A}(#1)} \newcommand{\sComp}[1]{{#1}^{c}} \newcommand{\sIn}[1]{{#1}^{i}} \newcommand{\sAd}[1]{{#1}^{a}} \newcommand{\sOp}[2]{{#1}^{#2}} \newcommand{\rt}{\mathbf{t}} \newcommand{\ReLU}{\mathrm{ReLU}} \newcommand{\rx}{\mathbf{x}} \newcommand{\rf}{\mathbf{f}} \newcommand{\cf}{\mathrm{f}} \newcommand{\rv}{\mathbf{v}} \newcommand{\ry}{\mathbf{y}} \newcommand{\rz}{\mathbf{z}} \newcommand{\card}{\mathrm{card}} \newcommand{\od}{\mathrm{d}} \newcommand{\deriv}[2]{\frac{\partial #1}{\partial #2}}

平面と粒子の衝突課題

Hiroki Shibata, Tokyo Metropolitan University
Creative Commons License version 4, CC BY, NC

はじめに

この課題は,template.pyのクラスR2, R2R2を拡張して取り組むものとする.R2は2次元ベクトル,R2R2は2次元ベクトル空間の上の線形変換(正方行列)のクラスである.基本的な加減乗がすでに実装してある.

実問題の題材として,平面と粒子の衝突課題を扱う.課題は6つからなる

課題1: 逆行列の実装

R2R2クラスにメンバ関数inverseを実装する.inverse は,呼び出しもとのオブジェクトが示す行列の逆行列をR2R2クラスのオブジェクトとして生成し,戻り値として返すように実装する.

逆行列関数の実装を追加したR2R2のオブジェクト(m とする)の各メンバ変数を,乱数で値を初期化し,その逆行列(mIとする)を求め,m*mI がほぼ単位行列となることを確認せよ.乱数による計算は,少なくとも32回繰り返さなくてはならない.具体的には,template.pyの「課題1 の試験コード」と書かれた部分が動くように関数 inverse を実装する.

m*mI の結果得られる行列を\(R\),この4つの要素を,\(R_{11}, R_{12}, R_{21}, R_{22}\)として,この順で要素の値を画面に出力した例を以下に示す.1行が一回の試行に対応する.e-16 などは,\(10^{-16}\)という意味であり,これがかかった値は浮動小数点数の誤差の範囲程度の小ささである.どの結果も,\(R\)がほぼ単位行列であることを示している.

(1.0,0.0:0.0,1.0)
(0.9999999999999998,0.0:-4.440892098500626e-16,0.9999999999999998)
(1.0,0.0:2.220446049250313e-16,1.0)
(1.0,0.0:0.0,0.9999999999999999)
(1.0,0.0:0.0,1.0)
(0.9999999999999999,0.0:0.0,0.9999999999999999)
(1.0,3.469446951953614e-18:-4.440892098500626e-16,1.0)
(1.0000000000000004,-2.220446049250313e-16:0.0,1.0)
(1.0,0.0:0.0,0.9999999999999999)
(1.0,0.0:0.0,1.0)
(0.9999999999999999,0.0:-5.551115123125783e-17,0.9999999999999999)
(1.0,0.0:-1.1102230246251565e-16,1.0)
(1.0,0.0:0.0,1.0)

課題2: 線分可視化

2点をR2のオブジェクト2つで定義し,その間を結ぶ線分から無作為に選んだ点の集合により,線分を可視化する.

2つのベクトル\(a, b\in\sR^2\)があるとする.\(a, b\)を結ぶ線分上の点\(v\)は次の式で表される. v = a t + b(1-t), t\in [0, 1] \(t\)を区間\([0, 1]\)の一様乱数により\(N\)生成し,この列を\(t_i, i..., N\)とする.このとき, v_i = a t_i + b(1-t_i), i=1,2,...,N のように,多数の線分上の点\(v_i\)が生成される.各\(i\)に対する\(v_i\)の座標を\(x_i, y_i\)と表す.\(x_i, y_i\)を点として描画し,直線を可視化する.

\(p_1 = (0.7,0.9), p_2 = (0.8,0.4), p_3 = (0.8, 0.7), p_4 = (0.2,0.9), p_5 = (0, 0.2), p_6 = (0.9,0.5)\) ,標本数を\(N=128\)とし,\(\overline{p_1p_2}, \overline{p_3p_4}, \overline{p_5p_6}\)を描画した例をに示す.ここで設定した以外の座標を持つ点を適宜選び,同様に可視化せよ.ただし,点の対は一つで良い.

3つの線分を描画した例

matplotlib を plt として,その表示領域について,以下,

plt.xlim(0, 1)
plt.ylim(0, 1)
を設定しておくことを推奨する.このように設定すると,点の分布範囲に関わらず,表示領域が\([0,1]\times[0,1]\)に固定される.今回の実験ではこのほうが結果を確認しやすい.見づらい場合は,設定値を変える.

課題3: 線分可視化の点集合を返す関数

課題2の手順を一部関数により実装し,汎用性を向上させる.具体的に,2つのR2のオブジェクトを受け取り,線分の点の第1座標のリストと,同第2座標のリストの2つのリストを戻り値として返す関数を lineBetween として作成し,その機能を試験する.関数の仕様は次のとおりである.

template.pyの「# 課題3」と書かれた部分のように使うことができ,返り値のlax, lay を指定し,グラフを描画し,動作を確認せよ.

課題4: 当たり判定

ある開始点\(v_0\)から\(d\)の方向へベクトルを伸ばしていき,\(p_1, p_2\)を結ぶ線分に交差する場合は,その交差点で伸長を停止させ,交差しない場合は,十分長く伸ばす,という手順を実装する.これにより,直線上に動く粒子が壁にぶつかるかどうか,またそのぶつかる場所を計算することが可能となる.

以上を考えるための方程式は次のとおりである.左辺が,\(v_0\)を起点とし\(d\)の方向へ伸びる直線(これを粒子の軌跡とする),右辺が\(p_1, p_2\)を通る直線(壁と呼ぶ)を表している.\(a_1, a_2\)は媒介変数であり,二次元ベクトルの方程式であるため,この2つを求められる可能性がある. v_0 + d a_1 = (p_2 - p_1)a_2 + p_1 上記を移行して,行列を用い表すと次のとおりになる. v_0 - p_1 = A \begin{pmatrix}a_1\\a_2\end{pmatrix} ここで,\(A\)は\(p_2 - p_1\)を第1列,\(-d\)を第2列に持つ行列である.この方程式は,\(p_2 - p_1\)と\(d\)が並行出ない場合に,またそのときに限り,次のように解くことができる. \begin{pmatrix}a_1\\a_2\end{pmatrix} = A^{-1}(v_0 - p_1) 上の式により求まる\(a_1, a_2\)を調べれば,粒子が壁にぶつかるかどうかがわかる.具体的に, a_1 \gt 0, 0\leq a_2 \leq 1 であれば,粒子は壁を貫通する.粒子と壁がちょうど重なる座標はもちろん\(v_0 + d a_1\)あるいは\((p_2 - p_1)a_2 + p_1\)である.

\(v_0 = (0.1,0.1), d=(0.5,0.5), p_1=(0.8,0), p_2=(0.8,0.9)\)とした場合の結果をへ示す.ちょうど衝突するところで直線の伸びが停止していることがわかる.このとき,\((a_1, a_2) \approx (0.896, 0.333)\)となっている.係数の値もこのように一緒に確認すると間違いが有った場合に原因を特定しやすい.

壁への衝突の例

今示しめした設定値以外で,色々試せ.

課題5: 影の計算

課題4の\(d\)の方向を,乱数で決定し,いろいろな方向へ粒子を射出する様子をシミュレートする.壁に当たるまで粒子の軌跡を伸ばし,その結果を今までの方法で可視化する.以下のような結果が得られるように,コードを作成せよ.壁の大きさ,位置,向き,粒子の数や方向は自由に変えて良い.

影の計算結果の例

ここで,軌跡上の点の間隔が一定になる工夫をしたほうが結果が綺麗に映る.の結果はこれを施したものである.このために,式の答えの\(a_1\)が得られたら,点の数\(N\)を,ある定数\(C\in\sR, N_0\)を用いて, N = \mathrm{int}(N_0 a_1/C ) と設定すると良い.\(d\)が単位ベクトルなら,\(C\)は線分\(\overline{v_0 (v_0 + da_1)}\)の最大の長さを意味する.\(N_0\)は最大の長さとなるときの,点の数である(最大の点数である).

壁を赤い線で描画するために matplotlib の scatter を呼び出す他,軌跡の点に対する matplotlib の scatter の呼び出しは,すべての点のリストを,一つのリストにまとめて一度にするとよい.

課題6: 反射

粒子の軌跡を壁で反射させる方法を述べる.プログラムが得意だが線形代数が苦手なものは特に取り組まれよ.壁の2点を\(p_1, p_2\), 法線ベクトルを\(h\),この壁にあたった場所を\(v_1\),粒子の射出点を\(v_0\),射出方向を\(d_0\)とする.\(v_1\)は\((d_0, v_0, p_2, p_1)\)の組から式により定めるものとする.ここで,\(h\)について, h = M\br{\frac{\pi}{2}} \frac{p_2 - p_1}{\nbr{p_2-p_1}} と定義する.\(M(\theta)\)は\(M(\theta)q\)としたとき,\(q\)を原点を中心として\(\theta\)半時計回りに回転させる,回転行列であるとする.つまり,法線ベクトルとは,方向ベクトル\(p_2 - p_1\)を90度左に回転させたものである.3次元だと外積など使うが,2次元ではこうである.

により\(h\)を求めたら,これを用いて,\(v_1\)から反射した後射出される方向\(d_1\)は,ベクトル\(a, b\)の内積を\(\exbr{a, b}\)とかくとして, d_1 = d_0 - 2\exbr{d_0, h} h となる.に図を示してある.軌跡の法線方向に向かう成分を打ち消すのが,\(\exbr{d_0, h} h\)である.法線方向への射影をとり,法線ベクトルにかけて射影の長さを持った法線ベクトル方向のベクトルを得ている.これを\(d_0\)から引くだけだと,法線方向成分が打ち消されて,\(d_1\)は壁の接線方向へ向かうようになる.もう一度同じベクトルを引くことで,反対向きに反射するようになる.だからでは2倍している.逆にも考えてみよう.\((d_1 + d_0)/2\)は,\(d_1, d_2\)の中線となるはずだがこれが今,\(h\)の方向と一致しているのだから,\(h\)は反射の対称線になっているとわかる.

結果をに示す.この反射は,実世界ではほぼ起こらない完全なものである.その軌跡の美しさに感動する.

粒子と壁の反射の例

課題7: 多段反射

粒子の現在の移動距離を\(L_t\),最大で許容する移動距離を\(L_*\)とする.現在の粒子の位置を\(v_t\), 移動方向を\(d_t\)とする.\(\nbr{d_0}=1\)であるとする.壁の列\(p_i, i=1,2,...,2L\)を用意する.多重反射させる場合は,進んだ距離を図るために方向ベクトルの大きさが1である方がわかりやすいため,\(\nbr{d_t}=1\)としている.このとき,で計算される\(d_1\)は\(\nbr{d_1}=1\)を満たすので,最初の\(t=0\)に関する\(d_0\)だけ\(\nbr{d_0}=1\)と設定しておけば,ここで述べる反射において\(\nbr{d_t}=1\)がどの\(t\)でも保たれる.

繰り返しの段階\(t\)において,線分\(\overline{p_{2m-1} p_{2m}}\)(壁)と粒子の衝突をすべての\(m=1,2,...,L\)について調べ,式を計算し,最小の\(a_1\)を与える壁を求める.もとの反射の条件に加え,\(a_1 \leq L\)ならば反射を実行する.\(a_1 > L\)ならば,\(a_1 =L\)として,\(v_{t+1} = v_t + d_t a_1\)により\(v_{t+1}\)を求めて,線分\(\overline{v_{t}v_{t+1}}\)を描画しこの粒子への処理を就床する.

反射させる,\(a_1 \leq L\)の場合の処理を述べる.壁を突き抜けないように,少し手前で反射させるために,\(\beta a_1, \beta\approx 0.99, \beta \lt 1\)のように少し小さく修正された衝突までの距離を用いて,衝突点を\(v_{t+1}\)として, v_{t+1} = v_t + d_t a_1\beta を計算する.\(d_{t+1}\)はと同様に, d_{t+1} = d_{t} - 2\exbr{d_t, h_*} h_* と求める.ただし,\(h_*\)は反射させかつ最小の\(a_1\)を与えた壁(最近傍壁)の法線ベクトルであり,に最近傍壁の情報を適用して得た\(h\)である.

残り長さを L_{t+1} = L_t - a_1 \beta と更新して,この節の処理を繰り返す.\(L_{t+1} \le 10^{-9}\)程度の小さな誤差の範囲に入ってきたら,この粒子に対する処理を終了する.例をへ示す.

64の粒子の多重反射の例

おわりに

なぜ乱数で点を配置するか,それは,等間隔に配置しようとすれば,難しいからだ.配置するべき長さが,その等間隔の幅で割り切れないかもしれない.乱数で確率的に与えれば,こういったことに悩む必要はない.点の密度の期待値は,線分の長さに比例する.