pianofisica

Mathematics & Physics, Maxima, a bit Python & Wolfram, and Arts

Python(SymPy)によるベクトル解析その1(ベクトルの内積と外積)

数学の具体的な計算にPython(SymPy)を使って、数学もPython(SymPy)も同時に学んでしまいましょう。今回はPython(SymPy)を使ってベクトル解析をしてみたいと思います。具体的にはベクトルの内積外積を定義し、その性質やいくつかの公式をみてみます。

Python(SymPy)のごく基本的な使い方については以下の記事を参照してください:

pianofisica.hatenablog.com



ベクトルの入力

ベクトル

 \displaystyle{\vec{A}=\left(\begin{array}{c} a_1 \\ a_2 \\ a_3 \end{array}\right) \qquad \vec{B}=\left(\begin{array}{c} b_1 \\ b_2 \\ b_3 \end{array}\right) \qquad \vec{C}=\left(\begin{array}{c} c_1 \\ c_2 \\ c_3 \end{array}\right)}

はSymPyでは(ひとつの方法として)以下のように入力します。

import sympy as sp
sp.var(' a1:4 b1:4 c1:4 ')
A = sp.Matrix([[a1],[a2],[a3]])
B = sp.Matrix([[b1],[b2],[b3]])
C = sp.Matrix([[c1],[c2],[c3]])

ベクトルの和、スカラー

ベクトル同士の和は、ベクトルの各成分の和として定義されます。つまり

 \displaystyle{\vec{A}+\vec{B}=\left(\begin{array}{c} a_1+b_1 \\ a_2+b_2 \\ a_3+b_3 \end{array}\right)}

とします。したがってベクトルの和は交換法則を満たします:

 \displaystyle{\vec{A}+\vec{B}=\vec{B}+\vec{A}}

またベクトルの和は結合法則も満たします:

 \displaystyle{(\vec{A}+\vec{B})+\vec{C}=\vec{A}+(\vec{B}+\vec{C})}

ベクトルの和の定義により、すべての成分がゼロであるゼロベクトル

 \displaystyle{\vec{0}=\left(\begin{array}{c} 0 \\ 0 \\ 0 \end{array}\right)}

は、ベクトルの和の演算に関する単位元であることがわかります。つまり

 \displaystyle{\vec{A}+\vec{0}=\vec{0}+\vec{A}=\vec{A}}

です。ベクトルの和に対するベクトル  \vec{A} の逆元  -\vec{A} は、ベクトルのすべての成分についてマイナスをとったものとして定まります。

 \displaystyle{\vec{A}+(-\vec{A})=\vec{0}}

いまスカラー k に対して、ベクトル  \vec{A} k 倍を、ベクトルの各成分の  k 倍として定義します。

 \displaystyle{k\vec{A}=\left(\begin{array}{c} ka_1 \\ ka_2 \\ ka_3 \end{array}\right)}

これは、ベクトルの和に対するベクトル  \vec{A} の逆元を  -\vec{A} と表したことと整合します。以上から、ベクトルの演算について、 \vec{A}+(-\vec{B}) \vec{A}-\vec{B} と表しても問題がないことがわかります。

さて、SymPyにおけるベクトルの和とスカラー倍ですが、これらは以下のようにして求められます。

A+B
sp.var(' k ')
k*A
A+(-1*A)

基底ベクトルの入力

以上のベクトルの和とスカラー倍の考察から、ベクトルの一般の元は以下の基底ベクトルの線形結合として表せることがわかります。つまり、基底ベクトルの組を

 \displaystyle{\vec{e}_1=\left(\begin{array}{c} 1 \\ 0 \\ 0 \end{array}\right) \qquad \vec{e}_2=\left(\begin{array}{c} 0 \\ 1 \\ 0 \end{array}\right) \qquad \vec{e}_3=\left(\begin{array}{c} 0 \\ 0 \\ 1 \end{array}\right)}

とするとき、ベクトル  \vec{A}

 \displaystyle{\vec{A}=a_1\vec{e}_1+a_2\vec{e}_2+a_3\vec{e}_3}

として表すことができます。

E1 = sp.Matrix([[1],[0],[0]])
E2 = sp.Matrix([[0],[1],[0]])
E3 = sp.Matrix([[0],[0],[1]])
A == a1*E1+a2*E2+a3*E3


 

ベクトルの内積とノルム

ベクトル  \vec{A},\ \vec{B}内積  \vec{A}\cdot\vec{B}

 \displaystyle{\vec{A}\cdot\vec{B}= a_1b_1+a_2b_2+a_3b_3}

によって定めます。とくにベクトルのそれ自体との内積は(ベクトルの各成分が実数であるとき)

 \displaystyle{\vec{A}\cdot\vec{A}= (a_1)^2+(a_2)^2+(a_3)^2}

により非負値です。そこでこの右辺の2乗根をとることができますが、その正の2乗根をベクトルのノルムといい

 \displaystyle{|A|=\sqrt{\vec{A}\cdot\vec{A}}}

で表すことにします。ベクトルの内積をSymPyで計算するには(あとの都合上、関数として定義しておきます)

def Dot(A, B): return sp.expand((A.transpose()*B)[0,0])
Dot(A,B)

とします。ノルムの計算は

def Nor(A): return sp.sqrt(Dot(A,A))
Nor(A)

とします。

内積の定義から、内積は2つのベクトルについて対称的な演算です。つまり

 \displaystyle{\vec{A}\cdot\vec{B}=\vec{B}\cdot\vec{A}}

が成り立ちます。また内積は(係数が実数であるとき)線形な演算であることもわかります。つまり

 \displaystyle{\vec{A}\cdot \left(k_1\vec{B}+k_2\vec{C}\right)=k_1\left(\vec{A}\cdot\vec{B}\right)+k_2\left(\vec{A}\cdot\vec{C}\right)}

が成り立ちます。これらの性質が成り立っているかSymPyで具体的に確認してみましょう。

Dot(A,B)==Dot(B,A)
sp.var(' k1 k2 ')
Dot(A,k1*B+k2*C)==sp.expand(k1*Dot(A,B)+k2*Dot(A,C))

2つのベクトルがなす角度

2つのベクトル  \vec{A},\ \vec{B} がなす角度  \varphi_{\vec{A}\vec{B}}

 \displaystyle{ \varphi_{\vec{A}\vec{B}}=\arccos\left(\frac{\vec{A}\cdot\vec{B}}{|\vec{A}|\,|\vec{B}|}\right) }

によって定義します。SymPyで計算するには

def Ang(A, B): return sp.acos(Dot(A,B)/(Nor(A)*Nor(B)))
Ang(A,B)

とします。定義から  \varphi_{\vec{A}\vec{B}}=\varphi_{\vec{B}\vec{A}} です。また

 \displaystyle{ \varphi_{\vec{A}\vec{A}}=\arccos\left(\frac{\vec{A}\cdot\vec{A}}{|\vec{A}|^2}\right)=\arccos(1)=0}

であることもわかります。

いくつか具体的にみてみましょう。

Ang(E1,E1)
Ang(E1,E2)
sp.var(' theta ', real=True, range=(0,sp.pi) )
sp.trigsimp(Ang(E1,sp.cos(theta)*E1+sp.sin(theta)*E2))

ベクトルの外積

2つのベクトル  \vec{A},\ \vec{B} から新たな第3のベクトルを得る演算のひとつとして、次で定義されるベクトルの外積があります。

 \displaystyle{\vec{A}\times\vec{B}= \left(\begin{array}{c} a_2b_3-a_3b_2 \\ a_3b_1-a_1b_3 \\ a_1b_2-a_2b_1 \end{array}\right) }

この演算の大きな特徴のひとつとして、積が交換法則を満たさないことがあげられます。実際

def Crs(A, B): return sp.expand(sp.Matrix([
[A[1,0]*B[2,0]-A[2,0]*B[1,0]],
[A[2,0]*B[0,0]-A[0,0]*B[2,0]],
[A[0,0]*B[1,0]-A[1,0]*B[0,0]]]))

によってSymPyでベクトルの外積を定義して、交換法則の成否を真偽判定してみると

Crs(A,B)==Crs(B,A)

より  \vec{A}\times\vec{B}\neq\vec{B}\times\vec{A} です。より詳細には、ベクトルの外積は反対称な積(順序を入れ替えると符号が変わる積)です。

 \displaystyle{\vec{A}\times\vec{B}=-\vec{B}\times\vec{A} }

実際、SymPyで確かめてみましょう。

Crs(A,B)==-Crs(B,A)

とくにベクトル  \vec{A} 自体の外積はゼロです。

 \displaystyle{\vec{A}\times\vec{A}=\vec{0} }

ベクトルの外積は、交換法則だけでなく、結合法則も満たしません。つまり

 \displaystyle{\vec{A}\times\left(\vec{B}\times\vec{C}\right)\neq\left(\vec{A}\times\vec{B}\right)\times\vec{C} }

です。実際

sp.expand(Crs(A,Crs(B,C)))==sp.expand(Crs(Crs(A,B),C))

です。したがってベクトルの外積に関する計算では、通常の数のかけ算のように、勝手なところから計算していいわけではありません。しかしながら、上にあげた積の反対称性と、次にあげる線形性は、積の性質を特徴付けています。

 \displaystyle{\vec{A}\times\left(k_1\vec{B}+k_2\vec{C}\right)=k_1\left(\vec{A}\times\vec{B}\right)+k_2\left(\vec{A}\times\vec{C}\right) }

これが成り立つことをSymPyで確かめてみましょう。

Crs(A,k1*B+k2*C)==sp.expand(k1*Crs(A,B)+k2*Crs(A,C))

ヤコビ恒等式

ベクトルの外積を特徴づけるもう一つの性質は、次のヤコビ恒等式(Jacobi identity)です。

 \displaystyle{\vec{A}\times\left(\vec{B}\times\vec{C}\right)+\vec{B}\times\left(\vec{C}\times\vec{A}\right)+\vec{C}\times\left(\vec{A}\times\vec{B}\right)=0}

積の反対称性を使ってこの等式を書き直すと

 \displaystyle{\vec{A}\times\left(\vec{B}\times\vec{C}\right)-\left(\vec{A}\times\vec{B}\right)\times\vec{C}=\left(\vec{C}\times\vec{A}\right)\times\vec{B}}

です。上のほうで、ベクトルの外積結合法則を満たしていないことをみましたが、この式は、ベクトルの外積結合法則の壊れ具合を定量的に示している式だと理解することができます。

ヤコビ恒等式が成り立っていることをSymPyで確かめてみましょう。

Crs(A,Crs(B,C))-Crs(Crs(A,B),C)==Crs(Crs(C,A),B)

実は、一般に対象に定められた演算が、1:反対称性、2:線形性、3:ヤコビ恒等式、という3つの条件を満たすとき、この演算はリー代数(Lie algebra)であるといいます。つまり、3次元の実ベクトルは、ベクトルの外積についてリー代数になっているのです。

単位ベクトルに関する公式、クロネッカーのデルタ、レヴィ・チヴィタ記号

さて、ベクトルの内積外積の単位ベクトル  \vec{e}_i\ (i=1,2,3) への作用をみてみましょう。というのは、内積外積も、線形性を備えた演算だったので、基底ベクトルへの作用さえわかれば一般のベクトルはそれらの線形結合として表されているすぎないためです。単位ベクトル  \vec{e}_i\ (i=1,2,3) について以下が成り立ちます。まずは内積

 \displaystyle{\begin{aligned} &\vec{e}_1\cdot\vec{e}_1=1\qquad\vec{e}_1\cdot\vec{e}_2=0\qquad\vec{e}_1\cdot\vec{e}_3=0\\ &\vec{e}_2\cdot\vec{e}_1=0\qquad\vec{e}_2\cdot\vec{e}_2=1\qquad\vec{e}_2\cdot\vec{e}_3=0\\ &\vec{e}_3\cdot\vec{e}_1=0\qquad\vec{e}_3\cdot\vec{e}_2=0\qquad\vec{e}_3\cdot\vec{e}_3=1\end{aligned}}

次に外積

 \displaystyle{\begin{aligned} &\vec{e}_1\times\vec{e}_1=\vec{0}\qquad\vec{e}_1\times\vec{e}_2=\vec{e}_3\qquad\vec{e}_1\times\vec{e}_3=-\vec{e}_2\\ &\vec{e}_2\times\vec{e}_1=-\vec{e}_3\qquad\vec{e}_2\times\vec{e}_2=\vec{0}\qquad\vec{e}_2\times\vec{e}_3=\vec{e}_1\\ &\vec{e}_3\times\vec{e}_1=\vec{e}_2\qquad\vec{e}_3\times\vec{e}_2=-\vec{e}_1\qquad\vec{e}_3\times\vec{e}_3=\vec{0}\end{aligned}}

規則性が見えるでしょうか?これらの関係は、より簡潔に次のように表すことができます。

 \displaystyle{\vec{e}_i\cdot\vec{e}_j=\delta_{i,j}\qquad\quad\vec{e}_i\times\vec{e}_j=\sum_{k=1}^3\epsilon_{ijk}\vec{e}_k}

ここで  \delta_{i,j}クロネッカーのデルタ(Kronecker delta)で、添字  i,\ j が等しいときには1を、異なるときにはゼロを与える関数です。また  \epsilon_{ijk} はレヴィ・チヴィタ記号(Levi-Civita symbol)といい、添字  i,\ j,\ k が1、2、3の偶置換の場合には1を、奇置換の場合にはマイナス1を、それ以外の場合にはゼロを与える関数です。これらの関数はSymPyでも組み込み関数として用意されています。

クロネッカーのデルタ

sp.KroneckerDelta(1,1)
sp.KroneckerDelta(1,2)

レヴィ・チヴィタ記号

sp.LeviCivita(1,2,3)

または

sp.Eijk(1,3,2)

この記事では3次元のベクトル空間を考えているので、添字として1、2、3しかとらない場合を考えますが、これらの組み込み関数はより多くの添字の場合にも対応しています。例えば

sp.LeviCivita(1,2,3,4,5,6,7,8)

いくつかの有用な公式

ベクトルの内積外積について以下の公式が成り立ちます。

 \displaystyle{\begin{aligned} &\vec{A}\cdot(\vec{B}\times\vec{C})=\vec{B}\cdot(\vec{C}\times\vec{A})\\ &\vec{A}\times(\vec{B}\times\vec{C})=\vec{B}(\vec{A}\cdot\vec{C})-\vec{C}(\vec{A}\cdot\vec{B})\\ &(\vec{A}\times\vec{B})\cdot(\vec{C}\times\vec{D})=(\vec{A}\cdot\vec{C})\,(\vec{B}\cdot\vec{D})-(\vec{A}\cdot\vec{D})\,(\vec{B}\cdot\vec{C})\\&(\vec{A}\times\vec{B})\times(\vec{C}\times\vec{D})=(\vec{A}\cdot(\vec{C}\times\vec{D}))\vec{B}-(\vec{B}\cdot(\vec{C}\times\vec{D}))\vec{A}\end{aligned}}

SymPyでこれらの公式を確かめてみましょう。

import sympy as sp
sp.var(' a1:4 b1:4 c1:4  d1:4 ')
A = sp.Matrix([[a1],[a2],[a3]])
B = sp.Matrix([[b1],[b2],[b3]])
C = sp.Matrix([[c1],[c2],[c3]])
D = sp.Matrix([[d1],[d2],[d3]])
def Dot(A, B): return sp.expand((A.transpose()*B)[0,0])
def Crs(A, B): return sp.expand(sp.Matrix([
[A[1,0]*B[2,0]-A[2,0]*B[1,0]],
[A[2,0]*B[0,0]-A[0,0]*B[2,0]],
[A[0,0]*B[1,0]-A[1,0]*B[0,0]]]))
Dot(A,Crs(B,C))==Dot(B,Crs(C,A))
Crs(A,Crs(B,C))==sp.expand(Dot(A,C)*B-C*Dot(A,B))
Dot(Crs(A,B),Crs(C,D))==sp.expand(Dot(A,C)*Dot(B,D)-Dot(A,D)*Dot(B,C))
Crs(Crs(A,B),Crs(C,D))==sp.expand(Dot(A,Crs(C,D))*B-Dot(B,Crs(C,D))*A)

 

キーワード:ベクトル解析、SymPy、Python