#!/usr/bin/env python # coding: utf-8 # # Hiroshi TAKEMOTO # (take.pwave@gmail.com) #

ベクトルと行列

#

# ベクトルと行列の計算は、幅広い分野で使われています。 # Sageを使ってベクトルと行列の特徴的な計算方法について見ていきましょう。 #

# #

ベクトルの生成

#

# sageでベクトルを生成する関数は、vector関数です。vector関数の使い方を以下に示します。 # 

# vector(値のリスト)
# または、
# vector(環, 値のリスト)
#
# #

#

# ベクトルの生成例として、$v=(2, 1, 3), w=(1, 1,-4)$の例を以下に示します。環については後ほど紹介します。 #

# # In[33]: # jupyter用のdisplayメソッド from IPython.display import display, Latex, HTML, Math, JSON # In[1]: # ベクトルと行列 v = vector([2,1,3]); v # In[2]: w = vector([1,1,-4]); w # In[3]: # ベクトルの和 v+w # #

ベクトルの計算

#

# Sageではベクトルの和と差は通常の変数と同じように+, -記号で表すこととができます。 #

#

# ベクトルのスカラー倍は、ベクトルの各要素にスカラを掛けた形になります。 #

#

# ベクトルの絶対値は、ベクトル用各要素の自乗の和を平方根となり、ベクトルの距離は２つのベクトルの差の絶対値となります。 #

#

# このようにSageではベクトルに対して、変数と同じように計算ができる点が特徴です。 #

# # In[4]: var('a1 a2 a3 b1 b2 b3 c') A = vector([a1, a2, a3]) B = vector([b1, b2, b3]) show (A+B) # ベクトルの和 show (c*A) # ベクトルのスカラー倍 show (abs(A)) # ベクトルの絶対値 show (abs(A-B)) # ベクトルの距離 # #

ベクトルの内積

#

# ベクトルには、内積と外積という特徴的な２つの演算があります。 #

#

# 内積は、２つのベクトルの各要素の積の和として、以下のように定義されます。 # $$ # \mathbf{v}\cdot\mathbf{w} = \Sigma_{i=1}^N v_i w_i # $$ #

#

# 先に生成したベクトルA, Bに対して内積を取ると、定義どおりの結果が得られます。 #

# # In[6]: AdB = A.dot_product(B) show(AdB) # #

# 内積の大きな特徴に、２つのベクトルのなす角度$\theta$との間に以下の関係があります。 # $$ # \mathbf{v}\cdot\mathbf{w} = |\mathbf{v}||\mathbf{w}| cos \theta # $$ # この性質を使って、ベクトルの類似を計算するのに、$cos(\theta)$がよく使われます。 #

#

# また、２つのベクトルが直行する場合には、$cos \theta = 0$から、 # ベクトルの内積はベクトルの直交判定にもよく利用されます。 #

#

# 以下に、$v=(2, 1, 3), w=(1, 1,-4)$の内積の結果とベクトルv, wのなす角度$\theta$を # 求めています。 #

# # In[7]: v.dot_product(w) # In[8]: deg = lambda x : N(x * 180/pi) # vとwからcos(th)を計算 cos_th = v.dot_product(w)/(abs(v)*abs(w)) th = arccos(cos_th) print deg(th) # #

# ３次元プロットで、ベクトルv,wをプロットしたのが以下の図です。 # ２つのベクトルのなす角度が約120度であることが、見て取れます。 #

# # In[10]: ZERO = vector([0, 0, 0]) v_line = line([ZERO, v], rgbcolor='blue') w_line = line([ZERO, w], rgbcolor='green') vw_line = line([ZERO, v.cross_product(w)], rgbcolor='red') (v_line+w_line).show(aspect_ratio=1) # #

ベクトルの外積

#

# 同様にベクトルの外積$\mathbf{v}\times\mathbf{w}$は、cross_product関数で計算します。 #

#

# ３次元ベクトルAとBの外積の結果は、以下のようなります。 #

# # In[11]: AcB = A.cross_product(B) show(AcB) # #

# ちょっと覚えにくいので、以下の性質を使って覚えると良いでしょう。 # $$ # \mathbf{A}\times\mathbf{B} = # \left| \begin{array}{ccc} # i & j & k \\ # a_{1} & a_{2} & a_{3} \\ # b_{1} & b_{2} & b_{3} # \end{array} \right| # $$ #

#

# x, y, z方向の単位ベクトルU=(i, j, k)を定義し、U, A, Bからなる行列mを作り、 # そのdetを求める、結果をそれぞれ、$i, j, k$でまとめると、 # $$ # (a_2 b_3 - a_3 b_2) i + (-a_1 b_3 + a_3 b_1) j + (a_1 b_2 - a_2 b_1) k # $$ # と先の外積の結果と対応することが分かります。 #

# # In[12]: var('i j k') U = vector([i, j, k]) m = matrix([U, A, B]); show(m) dm = det(m); show(expand(dm)) # #

# ベクトルの外積の図形的特徴は、外積の方向はベクトルvからベクトルwにねじを回して進む方向となり、 # その大きさはベクトルvとベクトルwの作る平行四辺形の面積となります。 # $$ # |\mathbf{v}\times\mathbf{w}| = |\mathbf{v}||\mathbf{w}|sin(\theta) # $$ #

#

# Sageを使ってベクトルv, wの外積を求め、その値が $|v||w| cos (\theta)$と一致することを見てみましょう。 #

# # In[13]: VcW = v.cross_product(w) # thは内積で求めた結果を利用 print VcW, N(abs(VcW)), N(abs(v)*abs(w)*sin(th)) # #

# ３次元プロットで、ベクトルv,wと外積の結果をプロットしたのが以下の図です。 # 図を少し回転すると、ベクトルvからベクトルwに回したときのねじの進行方向に外積ベクトルがのびている # ことが確認できます。 #

# # # In[14]: vw_line = line([ZERO, v.cross_product(w)], rgbcolor='red') # 図を回転しないと確認しずらいが、vからwにねじを回した方向になっている (v_line+w_line+vw_line).show(aspect_ratio=1) # #

行列の生成

#

# 行列は、matrix関数で生成します。 #

# 
# 	matrix(行列の要素のリスト)
# 	または
# 	matrix(環, 行列の要素のリスト)
#
#

# 生成された行列の要素は、配列と同じようにカギ括弧[]に要素のインデックスを指定することで、 # 参照できます。行を取得する場合には、行のインデックスを指定し、列を取得するにはcolumnメソッド # 列のインデックスを指定することで所望の情報を得ることができます。 #

# # In[15]: M = matrix([[1,2,3],[3,2,1],[1,2,1]]) show(M) # In[16]: print M[1] # 2行目（インデックスでは1）を取得 print M[0, 2] # 1行目3列の要素を取得 print M.column(1) # 2列目を取得 # #

# 行列とベクトルの積は、*演算子を使って行い、その結果としてベクトルが返されます。 #

# # In[18]: Mw = M*w; print M*w type(Mw) # #

行列の基本演算

#

# 行列の基本計算をSageを使ってみてみましょう。行列AとBの和、スカラー倍、積の結果を # 以下に示します。行列の積は、順序によって結果が異なることに注意してください。 #

# # In[19]: var('a11 a12 a21 a22 b11 b12 b21 b22') A = matrix([[a11, a12], [a21, a22]]) B = matrix([[b11, b12], [b21, b22]]) show(A) show(B) show(A+B) # 行列の和 show(c*A) # 行列のスカラー倍 show(A*B) # 行列の積（AB順） show(B*A) # 行列の積（BA順） # #

単位行列

# 単位行列は、identity_matrixで生成します。identity_matrixの引数には、行列の次数を指定します。 # # In[21]: Im = identity_matrix(3) show(Im) # #

対角行列

# 対角行列は、diagonal_matrixで生成します。 # 
# diagonal_matrix(対角要素のリスト)
#
# # # In[22]: D = diagonal_matrix([1, 2, 3]) show(D) # #

行列の解

#

# 行列MにベクトルXを掛けて、ベクトルYとなる場合、行列MとベクトルYからXを求めるメソッドがsolve_rightです。 # $$ # \mathbf{M}\mathbf{X} = \mathbf{Y} # $$ #

#

# solve_rightの例を以下に示します。 #

# # In[23]: M = matrix([[1,2,3],[3,2,1],[1,2,1]]); show(M) Y = vector([0,-4,-1]); show(Y) # In[24]: # solve_rightを使って求める X = M.solve_right(Y); view(X) # octaveの左除算オペレータ\を使って求める X = M \ Y show(X) # # $\mathbf{M}\mathbf{X}$を計算すると、ベクトルY(0, -4, -1)となることから解が正しいことが確認できます。 # # In[26]: M*X # #

連立方程式と行列の解の関係

#

# 先の行列式は、以下のような連立方程式と一致します。 # $$ # \left\{\begin{array}{rrr} # x_1 + 2 x_2 + 3 x_3 & = & 0 \\ # 3 x_1 + 2 x_2 + x_3 & = & -4 \\ # x_1 + 2 x_2 + x_3 & = & -1 # \end{array}\right. # $$ #

#

# 上記の連立方程式をsolve関数で計算結果とXの値が同じになることをSageで確かめてみましょう。 #

# # # In[27]: (x1, x2, x3) = var('x1,x2, x3') eq = [ [ x1 + 2*x2 + 3*x3 == 0], [3*x1 + 2*x2 + x3 == -4], [x1 + 2*x2 + x3 == -1]] sol = solve(eq, [x1, x2, x3]); show(sol) # #

転置行列

#

# 転置行列は、transpose関数で取得できます。 #

#

# 転置行列の性質は、 #

# 最後の性質は、行列の掛け合わせる順序を入れ替えるときに便利です。 #

# # In[28]: At = A.transpose(); show(At) # 転置行列 Bt = B.transpose(); show(Bt) # In[38]: print("(A^T)^T = A") show(At.transpose()) print("(A+B)^T = A^T + B^T") show(At + Bt) show((A+B).transpose()) # In[39]: print("(AB)^T = B^T A^T") show((A*B).transpose()) show(At*Bt) # #

行列式

#

# 行列式detは、以下のように定義されます。 # $$ # det(\mathbf{A}) = # \left| \begin{array}{cc} # a_{11} & a_{12} \\ # a_{21} & a_{21} # \end{array} \right| # $$ #

#

# 行列式は、逆行列の計算に使われるので、よく逆行列が存在するか否かの判別に使用されます。 #

#

# Sageでは、行列式はdet関数を使って計算されます。 #

# # # In[40]: show(A.det()) # In[41]: Mdet = M.det(); print Mdet # # #

逆行列

#

# 逆行列は、inverse関数で取得できます。 #

#

# 逆行列の性質は、以下の通りです。 #

#

# # # In[43]: Minv = M.inverse() show(Minv) # In[44]: print("$(A^{-1})^{-1} = A$") show(Minv.inverse()) # In[46]: print("(A^{T})^{-1} = (A^{-1})^{T}") show(M.transpose().inverse()) show(M.inverse().transpose()) # #

環

#

# 環は、数の集合を表します。 # よく使われる環を以下に示します。 #