pianofisica

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

Python(SymPy)によるベクトル解析その2(勾配・発散・回転)

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

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

pianofisica.hatenablog.com


ベクトル値関数の入力

3次元空間座標(3次元デカルト座標 (x,y,z) の関数  f=f(x,y,z) はSymPyでは以下のように入力します。

import sympy as sp
sp.var(' x y z ')
from sympy.core.function import Function
f = Function('f')(x,y,z)

したがってベクトルの各成分  v_i が関数であるようなベクトル値関数

 \displaystyle{\vec{V}=\left(\begin{array}{c} v_1 \\ v_2 \\ v_3 \end{array}\right)}

はSymPyでは以下のようにして入力します。

v1 = Function('v1')(x,y,z)
v2 = Function('v2')(x,y,z)
v3 = Function('v3')(x,y,z)
V = sp.Matrix([[v1],[v2],[v3]])

勾配(Gradient)

3次元の空間座標  (x,y,z) の関数  f に対して、ベクトル値関数を以下で与える演算を勾配(Gradient)といい  {\rm grad}\ {f} などと表します。

 \displaystyle{{\rm grad}\ {f}=\left(\begin{array}{c} \frac{\partial f}{\partial x} \\ \frac{\partial f}{\partial y} \\ \frac{\partial f}{\partial z} \end{array}\right)}

発散(Divergence)

3次元の空間座標  (x,y,z) のベクトル値関数  \vec{V} に対して、関数を以下で与える演算を発散(Divergence)といい  {\rm div}\ \vec{V} などと表します。

 \displaystyle{{\rm div}\ \vec{V}=\frac{\partial v_1}{\partial x}+\frac{\partial v_2}{\partial y}+\frac{\partial v_3}{\partial z} }

回転(Rotation/Curl

3次元の空間座標  (x,y,z) のベクトル値関数  \vec{V} に対して、新たなベクトル値関数を以下で与える演算を回転(Rotation/Curl)といい  {\rm rot}\ \vec{V} {\rm curl}\ \vec{V} などと表します。

 \displaystyle{{\rm rot}\ \vec{V}=\left(\begin{array}{c} \frac{\partial v_3}{\partial y}-\frac{\partial v_2}{\partial z} \\ \frac{\partial v_1}{\partial z}-\frac{\partial v_3}{\partial x} \\ \frac{\partial v_2}{\partial x}-\frac{\partial v_1}{\partial y} \end{array}\right)}

これらの演算は、(形式的な)ベクトル  \nabla

\nabla=\displaystyle{\left(\begin{array}{c} \frac{\partial}{\partial x} \\ \frac{\partial}{\partial y} \\ \frac{\partial}{\partial z} \end{array}\right)}

によって導入したとき、それぞれ、勾配:関数に作用、発散:ベクトル値関数に内積をとって作用、回転:ベクトル値関数に外積をとって作用したものであるとみなせます。つまり

 \displaystyle{\begin{aligned}&{\rm grad}\ {f}=\nabla\,f\\&{\rm div}\ \vec{V}=\nabla\cdot\vec{V}\\&{\rm rot}\ \vec{V}=\nabla\times \vec{V}\end{aligned}}

です。したがって前の記事

pianofisica.hatenablog.com

でみたようなベクトルの演算に関するさまざまな公式が成り立ちますが、  \nabla は厳密な意味でのベクトルではなくあくまでも微分作用素であり、関数やベクトルに作用してはじめて意味がある対象であることを忘れてはいけません。たとえば、ベクトルの外積について

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

が成り立ちますが、もちろんここで

 \displaystyle{{\rm rot}\ \vec{V}=\nabla\times \vec{V}\neq-\vec{V}\times\nabla}

です。

さて以上の演算をSymPyで実装してみます。勾配、発散、回転はそれぞれ以下のようにして定義できます。

def grad(f): return sp.Matrix([
[sp.diff(f,x)],[sp.diff(f,y)],[sp.diff(f,z)]])
grad(f)
def div(V): return sp.diff(V[0,0],x)+sp.diff(V[1,0],y)+sp.diff(V[2,0],z)
div(V)
def rot(V): return sp.Matrix([
[sp.diff(V[2,0],y)-sp.diff(V[1,0],z)],
[sp.diff(V[0,0],z)-sp.diff(V[2,0],x)],
[sp.diff(V[1,0],x)-sp.diff(V[0,0],y)]])
rot(V)

ラプラシアン

以下で定義される2階微分演算子ラプラシアン(Laplacian)といいます。

 \displaystyle{\triangle=\frac{\partial^2}{\partial x^2}+\frac{\partial^2}{\partial y^2}+\frac{\partial^2}{\partial z^2} }

関数  f に対しては

 \displaystyle{\triangle\,f=\frac{\partial^2 f}{\partial x^2}+\frac{\partial^2 f}{\partial y^2}+\frac{\partial^2 f}{\partial z^2} }

ベクトル値関数  \vec{V} に対しては

 \displaystyle{\triangle\,\vec{V}=\left(\begin{array}{c} \triangle\,v_1 \\ \triangle\,v_2 \\ \triangle\,v_3 \end{array}\right) }

と定義します。

いくつかの公式

勾配、発散、回転について以下の公式が成り立ちます。

 \displaystyle{\begin{aligned}&\triangle\,f={\rm div}({\rm grad}\,f)  \\ &{\rm rot}({\rm grad}\ f)=\vec{0}\\ &{\rm div}({\rm rot}\ \vec{V})=0\\ & {\rm grad}({\rm div}\ \vec{V})-{\rm rot}({\rm rot}\ \vec{V})=\triangle \vec{V}\\ &{\rm grad}\ (fg)=f\,{\rm grad}\,g+g\,{\rm grad}\,f\\ &{\rm grad}(\vec{V}\cdot\vec{W})=\vec{V}\times({\rm rot}\,\vec{W})+\vec{W}\times({\rm rot}\,\vec{V})+(\vec{V}\cdot\nabla)\vec{W}+(\vec{W}\cdot\nabla)\vec{V} \\&{\rm div}\, (f\vec{V})=f{\rm div}\ \vec{V}+({\rm grad}\,f)\cdot\vec{V} \\&{\rm div}\, (\vec{V}\times\vec{W})=\vec{W}\cdot({\rm rot}\ \vec{V})-\vec{V}\cdot({\rm rot}\ \vec{W}) \\ &{\rm rot}\ f\vec{V}=f\,{\rm rot}\,\vec{V}+({\rm grad}\ f)\times \vec{V}\\ &{\rm rot}\ (\vec{V}\times\vec{W})=\vec{V}({\rm div}\,\vec{W})-\vec{W}({\rm div}\,\vec{V})+(\vec{W}\cdot\nabla)\vec{V}-(\vec{V}\cdot\nabla)\vec{W}\end{aligned}}

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

sp.expand(div(grad(f)))
sp.expand(rot(grad(f)))
sp.expand(div(rot(V)))
sp.expand(grad(div(V))-rot(rot(V)))
g = Function('g')(x,y,z)
sp.expand(grad(f*g)-f*grad(g)-g*grad(f))
w1 = Function('w1')(x,y,z)
w2 = Function('w2')(x,y,z)
w3 = Function('w3')(x,y,z)
W = sp.Matrix([[w1],[w2],[w3]])
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]]]))
def DotNabla(A,f): return A[0,0]*sp.diff(f,x)+A[1,0]*sp.diff(f,y)+A[2,0]*sp.diff(f,z)
sp.expand(grad(Dot(V,W))-Crs(V,rot(W))-Crs(W,rot(V))-DotNabla(V,W)-DotNabla(W,V))
sp.expand(div(f*V)-f*div(V)-Dot(grad(f),V))
sp.expand(div(Crs(V,W))-Dot(W,rot(V))+Dot(V,rot(W)))
sp.expand(rot(f*V)-f*rot(V)-Crs(grad(f),V))
sp.expand(rot(Crs(V,W))-div(W)*V+div(V)*W-DotNabla(W,V)+DotNabla(V,W))

 

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