Python（SymPy）によるベクトル解析その１（ベクトルの内積と外積）

数学の具体的な計算に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}$

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

$\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)

とします。

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

$\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))

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

２つのベクトル $\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))

ベクトルの外積

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

$\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)

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

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

さて、ベクトルの内積、外積の単位ベクトル $\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$ が等しいときには１を、異なるときにはゼロを与える関数です。また $\epsilon_{ijk}$ はレヴィ・チヴィタ記号（Levi-Civita symbol）といい、添字 $i,\ j,\ k$ が１、２、３の偶置換の場合には１を、奇置換の場合にはマイナス１を、それ以外の場合にはゼロを与える関数です。これらの関数はSymPyでも組み込み関数として用意されています。

クロネッカーのデルタ

sp.KroneckerDelta(1,1)

sp.KroneckerDelta(1,2)

レヴィ・チヴィタ記号

sp.LeviCivita(1,2,3)

または

sp.Eijk(1,3,2)

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

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