Pythonの内積、外積、ノルムに関する計算

Pythonでベクトルの内積、外積、ノルム(大きさ)を計算してきましたが、これらの関連で興味深い計算をしてみます。

2つのベクトルのノルムと内積から角度を計算する

NumPyで2つのベクトルのなす角度を計算する

2つのベクトルn_aとn_bをnadarrayで定義し、それぞれのノルムと内積を計算します。一般に$\displaystyle \cos \theta = \frac{(n_a・n_b)}{\|n_a\|\|n_b\|}$の関係にあるので、内積をノルムの積で割ることで2つのベクトルのなす角のcos(コサイン)を計算することができます。

  1. #1 NumPyによる2つのベクトルの角度の計算
  2. import numpy as np
  3. n_a = np.array([1, 4, 1])
  4. n_b = np.array([2, 5, 5])
  5. normn_a = np.sqrt(sum(n_a**2))
  6. normn_b = np.sqrt(sum(n_b**2))
  7. dotn_ab = np.dot(n_a, n_b)
  8. cosn_ab = dotn_ab/(normn_a*normn_b)
  9. "n_aのノルム:{:.6f} n_bのノルム:{:.6f} 内積:{:.6f} cos_ab:{:.6f}".format(
  10. normn_a, normn_b, dotn_ab, cosn_ab)

n_aのノルム:4.242641  n_bのノルム:7.348469  内積:27.000000 cos_ab:0.866025

それぞれ小数点第6位まで表示するようにしましたが、いまひとつすっきりとしません。

SymPyで2つのベクトルのなす角度を計算する

SymPyで2つのベクトルの内積からなす角のcosを求める

SymPyモジュールを使うと、次の通り平方根についてルート付きで表示することができます。その前提として、SymPyとともに、IpythonのMathモジュールをインポートしておきます。このことにより、数式をわかりやすく表示することができます。

  1. # SymPyによるベクトルの角度の計算
  2. import sympy
  3. from IPython.display import Math
  4. s_a, s_b = sympy.symbols('s_a s_b')
  5. s_a = sympy.Matrix([1, 4, 1])
  6. s_b = sympy.Matrix([2, 5, 5])
  7. norm_a = s_a.norm(ord=2)
  8. norm_b = s_b.norm(ord=2)
  9. dot_ab = s_a.dot(s_b)
  10. cos_ab= dot_ab/(norm_a*norm_b)
  11. display(Math(r'n_aのノルム=%s' sympy.latex(norm_a)))
  12. display(Math(r'n_bのノルム=%s' % sympy.latex(norm_b)))
  13. display(Math(r'内積=%s' % sympy.latex(dot_ab)))
  14. display(Math(r'cos=%s' % sympy.latex(cos_ab)))

Sympyによる2つのベクトルの角度の計算
Sympyによる2つのベクトルの角度の計算

Matrixとしてベクトルを定義し、normメソッドでノルムを、dotメソッドで内積を計算して、内積を2つのノルムの積で割ることで、角度に対するcosの値をルート付きで表示することができます。このためIpython.displayのMath関数を使います。

SymPyでcosの値から角度を計算する

cosの値が求まったので、これを角度に変換するためarccos (アークコサイン) を計算します。このため、多くのコンピュータ言語と同様、SymPyでもarccosはacos関数を使って計算します。

  1. # cosの値から角度を計算(arccos)
  2. theta_ab=sympy.acos(cos_ab)
  3. display(Math(r'\theta=%s' % sympy.latex(theta_ab)))

$\displaystyle\theta=\frac{\pi}{6}$

SymPyでは角度はラジアン(radian:弧度法)で表示されますが、どうしもピンとこないので、角度(degree:度数法)による表示に変換します。

  1. # ラジアン(radian:弧度法)から角度(degree:度数法)への変換
  2. import math
  3. deg=math.degrees(sympy.acos(cos_ab))
  4. deg

30.000000000000004

mathモジュールのdegrees関数を使うのが最も簡単です。結果は次の通り、ぴったりとはいきませんが、ベクトルs_aとs_bがなす角度は30度であることがわかりました。

2つのベクトルがなす平行四辺形と三角形の面積の計算

2つのベクトルがなす平行四辺形と三角形の面積にはいろいろな求め方があるので、これらを振り返りながらPythonを使って計算します。

2つの辺の長さとこれをはさむ角度がわかっている場合

2つの辺の長さとこれをはさむ角度がわかっている場合、次の式で面積を計算します。

$平行四辺形の面積=aの辺の長さ×bの辺の長さ×sin(abをはさむ角度)$

という式になります。三角形の面積はこの1/2になります。

ここではすでに、aの辺の長さはnorm_a、bはnorm_b、2辺をなす角度はtheta_abで計算されているので、次のプログラムで計算することができます。

  1. # 四辺形、三角形の面積
  2. rec_1 = norm_a * norm_b * sympy.sin(theta_ab)
  3. tri_1 = rec_1/2
  4. display(Math(r'四辺形の面積=%s' % sympy.latex(rec_1)))
  5. display(Math(r'三角形の面積=%s' % sympy.latex(tri_1)))

結果は、平行四辺形の面積はrec_1、三角形の面積はtri_1となります。Math関数を使うと次の通り表示されます。

2辺の長さとそれをはさむ角度がわかっている場合
2辺の長さとそれをはさむ角度がわかっている場合

外積から計算する方法

2つのベクトルの各要素の値がわかっていれば、そこからなす平行四辺形の面積は、2つのベクトルの外積のノルムの値に一致します。ここではすでにs_aとs_bはsympy.matrixで定義されているので、ここから外積を求め、ノルムを取ることで平行四辺形と三角形の面積を計算します。

  1. # 外積から平行四辺形と三角形の面積を計算する
  2. cross_ab=s_a.cross(s_b)
  3. rec_2 = cross_ab.norm(ord=2)
  4. tri_2 = rec_2 /2
  5. display(Math(r'外積=%s' % sympy.latex(cross_ab)))
  6. display(Math(r'四辺形の面積=%s' % sympy.latex(rec_2)))
  7. display(Math(r'三角形の面積=%s' % sympy.latex(tri_2)))

結果は次の通りです。高校数学では外積を取り扱っていないようですが、このように応用範囲が広いのでとても残念な気がします。

外積から平行四辺形と三角形の面積を計算する
外積から平行四辺形と三角形の面積を計算する

2つのベクトルの終点を結ぶ直線を使い面積を計算する

ベクトルs_aの終点を始点として、ベクトルs_bの終点を終点とするベクトルs_abを取ります。このベクトルのノルムが三角形の3つめの辺の長さになります。

  1. # 2つのベクトルの終点を結ぶベクトルとそのノルム
  2. s_ab = s_b-s_a
  3. norm_ab=s_ab.norm(ord=2)
  4. display(Math(r'終点を結ぶベクトル=%s' % sympy.latex(s_ab)))
  5. display(Math(r'ベクトルのノルム=%s' % sympy.latex(norm_ab)))

ベクトルs_abとそのノルムは次の通りです。

2つのベクトルの終点を結ぶベクトル
2つのベクトルの終点を結ぶベクトル

また、2つの辺の長さがnorm_a、norm_bで、2つのベクトルをはさむ角度がcos_abで計算されているので、余弦定理を使うことにより3つ目の辺の長さを計算することができます。

  1. # 余弦定理による三角形の辺の計算
  2. sympy.sqrt(norm_a ** 2 + norm_b ** 2 -2*cos_ab*norm_a*norm_b)

答えは同じく\3\sqrt{2}\になります。

以上のことから、2つのベクトルとその終点を結ぶベクトルのノルムが求まりました。三角形の面積はヘロンの定理によって求めることができるので早速、計算してみます。

  1. # ヘロンの定理による三角形の面積の計算
  2. s=(norm_a+norm_b+norm_ab)/2
  3. tri_3=sympy.sqrt(s*(s-norm_a)*(s-norm_b)*(s-norm_ab))
  4. display(Math(r'三角形の面積=%s' % sympy.latex(tri_3)))
ヘロンの公式による三角形の面積の計算
ヘロンの公式による三角形の面積の計算

式が複雑になるので、計算結果も一見すると間違っているように見えます。この場合、SymPyではsimplify関数を使うと式を簡潔にすることができます。

  1. # 複雑な式を簡潔にする
  2. display(Math(r'簡潔にした表示=%s' % sympy.latex(sympy.simplify(tri_3))))

このことにより、前と同じ結果になります。

ヘロンの公式による計算結果を簡潔に表示
ヘロンの公式による計算結果を簡潔に表示

正弦定理や内積を使って三角形の面積を計算する

最後に、少し変わった方法を使って三角形の面積を計算します。

正弦定理を使う方法

2ベクトルがなす角度と、終点を結ぶベクトルのノルム(長さ)がわかると、三角形の外接円の半径を計算することができます。三角形の内角の正弦(サイン)とその対辺の長さの関係は辺の長さ/正弦が、外接円の直径(半径の2倍)になります。さらに、三角形の3辺の積を半径の4倍で割ることにより、三角形の面積を計算することができます。

  1. # 正弦定理を使った三角形の面積の計算
  2. r = norm_ab / (sympy.sin(theta_ab)*2)
  3. tri_4=norm_a * norm_b * norm_ab/(r*4)
  4. display(Math(r'三角形の面積=%s' % sympy.latex(tri_4)))

内積と使う方法

そもそも、2つのベクトルの要素がわかれば、外積が使えなくても内積を使って面積を計算することができます。三角形の面積は、2つの辺のノルムの2乗の合計から内積の2乗を差し引いた結果の平方根が面積になります。

  1. # 内積を使った三角形の面積の計算
  2. tri_5 = sympy.sqrt(norm_a **2 * norm_b **2 - dot_ab ** 2)/2b
  3. display(Math(r'三角形の面積=%s' % sympy.latex(tri_5)))

答えはいずれも$\frac{9\sqrt{3}}{2}$になります。

外積により6面体の体積を計算する

ベクトルs_a,s_bが張る平行六面体の体積は、次の式で計算することができます。

平行六面体の体積はs_aとs_bから計算した外積のベクトルとs_cの内積の絶対値となります。数式であらわすと次の計算をすることになります。

V=|(s_a×s_b)・s_c|

  1. # ベクトルの外積と内積から平行六面体の体積を計算する
  2. s_a = sympy.Matrix([1, 2, -1])
  3. s_b = sympy.Matrix([3, 1, 4])
  4. s_c = sympy.Matrix([1, 0, -1])
  5. cross_ab=s_a.cross(s_b)
  6. v_abc=cross_ab.dot(s_c)
  7. display(Math(r'外積=%s' % sympy.latex(cross_ab)))
  8. display(Math(r'体積=%s' % sympy.latex(v_abc)))

結果は同じです。