この課題は,template.pyのクラスR2, R2R2を拡張して取り組むものとする.R2は2次元ベクトル,R2R2は2次元ベクトル空間の上の線形変換(正方行列)のクラスである.基本的な加減乗がすでに実装してある.
実問題の題材として,平面と粒子の衝突課題を扱う.課題は6つからなる
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点をR2のオブジェクト2つで定義し,その間を結ぶ線分から無作為に選んだ点の集合により,線分を可視化する.
2つのベクトル\(a, b\in\sR^2\)があるとする.\(a, b\)を結ぶ線分上の点\(v\)は次の式で表される.
\(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}\)を描画した例をに示す.ここで設定した以外の座標を持つ点を適宜選び,同様に可視化せよ.ただし,点の対は一つで良い.
matplotlib を plt として,その表示領域について,以下,
plt.xlim(0, 1) plt.ylim(0, 1)を設定しておくことを推奨する.このように設定すると,点の分布範囲に関わらず,表示領域が\([0,1]\times[0,1]\)に固定される.今回の実験ではこのほうが結果を確認しやすい.見づらい場合は,設定値を変える.
課題2の手順を一部関数により実装し,汎用性を向上させる.具体的に,2つのR2のオブジェクトを受け取り,線分の点の第1座標のリストと,同第2座標のリストの2つのリストを戻り値として返す関数を lineBetween として作成し,その機能を試験する.関数の仕様は次のとおりである.
template.pyの「# 課題3」と書かれた部分のように使うことができ,返り値のlax, lay を指定し,グラフを描画し,動作を確認せよ.
ある開始点\(v_0\)から\(d\)の方向へベクトルを伸ばしていき,\(p_1, p_2\)を結ぶ線分に交差する場合は,その交差点で伸長を停止させ,交差しない場合は,十分長く伸ばす,という手順を実装する.これにより,直線上に動く粒子が壁にぶつかるかどうか,またそのぶつかる場所を計算することが可能となる.
以上を考えるための方程式は次のとおりである.左辺が,\(v_0\)を起点とし\(d\)の方向へ伸びる直線(これを粒子の軌跡とする),右辺が\(p_1, p_2\)を通る直線(壁と呼ぶ)を表している.\(a_1, a_2\)は媒介変数であり,二次元ベクトルの方程式であるため,この2つを求められる可能性がある.
\(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)\)となっている.係数の値もこのように一緒に確認すると間違いが有った場合に原因を特定しやすい.
今示しめした設定値以外で,色々試せ.
課題4の\(d\)の方向を,乱数で決定し,いろいろな方向へ粒子を射出する様子をシミュレートする.壁に当たるまで粒子の軌跡を伸ばし,その結果を今までの方法で可視化する.以下のような結果が得られるように,コードを作成せよ.壁の大きさ,位置,向き,粒子の数や方向は自由に変えて良い.
ここで,軌跡上の点の間隔が一定になる工夫をしたほうが結果が綺麗に映る.の結果はこれを施したものである.このために,式の答えの\(a_1\)が得られたら,点の数\(N\)を,ある定数\(C\in\sR, N_0\)を用いて,
壁を赤い線で描画するために matplotlib の scatter を呼び出す他,軌跡の点に対する matplotlib の scatter の呼び出しは,すべての点のリストを,一つのリストにまとめて一度にするとよい.
粒子の軌跡を壁で反射させる方法を述べる.プログラムが得意だが線形代数が苦手なものは特に取り組まれよ.壁の2点を\(p_1, p_2\), 法線ベクトルを\(h\),この壁にあたった場所を\(v_1\),粒子の射出点を\(v_0\),射出方向を\(d_0\)とする.\(v_1\)は\((d_0, v_0, p_2, p_1)\)の組から式により定めるものとする.ここで,\(h\)について,
式により\(h\)を求めたら,これを用いて,\(v_1\)から反射した後射出される方向\(d_1\)は,ベクトル\(a, b\)の内積を\(\exbr{a, b}\)とかくとして,
結果をに示す.この反射は,実世界ではほぼ起こらない完全なものである.その軌跡の美しさに感動する.
粒子の現在の移動距離を\(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}\)として,
残り長さを
なぜ乱数で点を配置するか,それは,等間隔に配置しようとすれば,難しいからだ.配置するべき長さが,その等間隔の幅で割り切れないかもしれない.乱数で確率的に与えれば,こういったことに悩む必要はない.点の密度の期待値は,線分の長さに比例する.