テンソルに親しもう
このワークシートはMath by Codeの一部です。
前回まで「慣性系」の相対性原理、つまり特殊な相対性理論を見てきたね。
ポイントは光速が絶対で一定であることだった。
相対性理論を計量まで深めるためには、時空のゆがみや変換を学ぶ必要があることに気づいたね。
そこで、時空を数式で表現するために、今回はテンソルについて学ぼう。
1。テンソルとは
＜テンソルは複数添え字（方向）をもつ数の入れ物＞
テンソルはデータが入れ物の型でベクトル・行列の1種ともいえます。
配列やリストが入れにすると作れます。
テンソルは自由度からみると、スカラー、ベクトルの次にくるもので、行列を広げたものです。
実際に、方向(軸）の数、階数（ランク）は
スカラーは０、ベクトルは１、テンソルは2以上です。
テンソルは、広い意味ではスカラー、ベクトルも含み、狭い意味では行列以上のものに限ります。
ベクトルと同様に、成分は座標系に連動して変わるけれどもテンソルそのものは変わらない。
・テンソルの型はshapeで表示できます。
・テンソルの軸はtransposeで,行列と同じように転置できます。
(例）
Python
#テンソルの型============
import numpy as np
scalar_x=np.array(2)
vector_x=np.array([1,2])
matrix_x=np.array([[1,2], [3,4]])
tensor_x=np.array([[[1,2], [3,4]], [[2,4], [6,8]]])
print(scalar_x.shape,",",vector_x.shape,",",matrix_x.shape,",",tensor_x.shape)
#================
[OUT}
() , (2,) , (2, 2) , (2, 2, 2)
(例）
1から24までの数の配列PをP.shape=(2,3,4)とすることで、
2×3×4のテンソルに整形しましょう。2段に区切り、3本に区切り、1本に4個ならべます。
順にz軸、y軸、x軸とすることで、ループをｚ、ｙを固定してｘだけ動かすと
ｘ軸にそって、１，２，３，４と数がならんだベクトルが1本できます。
それをｙ軸を動かしてから、ｘ軸にそって、５，６，７，８と数がならんだベクトルが1本できます。
ｚを動かすと、上の段にも数のベクトルを2本ならべられます。
・Q=P.transpose(0,2,1)とすることで、ｚ段数は変えずに行と列を転置できますね。
・Q=P.transpose(2,1,0)なら、ｙたて数を変えずに、ｚ段とｘよこの入れ替えができます。
Python
#テンソルの転置============
import matplotlib.pyplot as plt
import numpy as np
P=np.array(range(1,25))
P.shape=(2,3,4)
#段数をかえずに、たて・よこの入れ替え
Q=P.transpose(0,2,1)
Q
fig, ax = plt.subplots(1, 2, figsize=(10, 4), subplot_kw={'projection': '3d'})
def drawTensor(P,n):
print(P)
z_size, y_size, x_size = P.shape[0], P.shape[1],P.shape[2]
voxels = np.ones((x_size, y_size, z_size), dtype=bool)
ax[n].voxels(voxels, edgecolor='k', facecolors='green', alpha=0.1)
#テンソルの内容をその位置にかく
for pz in range(z_size):
for py in range(y_size):
for px in range(x_size):
contents = P[pz,py,px]
#print(f"{pz},{py},{px};contents")
ax[n].text(px, py, pz+0.5, f'{contents}', fontsize=15, verticalalignment="bottom")
ax[n].set_xlabel('x')
ax[n].set_ylabel('y')
ax[n].set_zlabel('z')
ax[n].set_aspect('equal')
drawTensor(P,0)
drawTensor(Q,1)
plt.show()
#====================================
[OUT]
[[[ 1 2 3 4]
[ 5 6 7 8]
[ 9 10 11 12]]
[[13 14 15 16]
[17 18 19 20]
[21 22 23 24]]]
[[[ 1 5 9]
[ 2 6 10]
[ 3 7 11]
[ 4 8 12]]
[[13 17 21]
[14 18 22]
[15 19 23]
[16 20 24]]]
２．テンソル計算をやってみよう
＜成分ごとの和差積商もある＞
テンソルはデータが入れ物の型でベクトル・行列の1種だから、
同じ位置のどうしのたし算、引き算、かけ算、わり算というものもできます。
A+B
A－B
A*B
A/B
これはベクトルの和と差、行列の和と差といっしょで、
同じ位置の成分どうしの演算なので、
テンソルの形状が同じものどうしの計算になり、
とくにテンソルらしい複雑さは何もありません。
なので、省略します。
＜テンソル積＞
テンソル積は、２つのテンソルで共通する方向（軸）を決めてその内積をとるものです。
共通軸を指定しないとき、
次元が q のテンソル A と
次元が r のテンソル B があるとき、
これらのテンソルの積C = A⊗Bは、 q + r 次元のテンソルになります。
・共通の方向（軸）の数を０とする(axes=0)と、２つのテンソルの各要素に対する直積になるので、
次元が和になります。１次元と１次元のテンソル積は1＋1＝２次元。
(例）
Python
[IN]
#============
from numpy import array
from numpy import tensordot
A = array([1,2])
B = array([3,4])
C = tensordot(A, B, axes=0)
print(C)
#=================
[OUT]
[[3 4]
[6 8]]
・共通の方向（軸）の数を１とする(axes=1)と、２つのテンソルの共通な軸で内積がおきるので、
その軸がつぶれます（縮約）。
軸数が１つかぶってそれぞれ、減ります。
２次元と２次元のテンソル積は行列の積と同じ2-1+2-1=2次元になります。
どの軸をつぶして内積をとるかをそれぞれのテンソルの添え字の番号で指定します。
(例）
Python
[IN]
#============
import numpy as np
# (2, 3) と (3, 4) の行列。共通の「3」の部分（軸1と軸0）で積をとる
A = np.array(range(1,7)) #１から６の数列
A.shape=(2,3)#2行3列の2軸
B = np.ones((3, 4))#１だけを3行4列の2軸
res = np.tensordot(A, B, axes=(1, 0))
resre = np.tensordot(B,A,axes=(0,1))
# 結果は (2, 4) の行列。np.dot(A, B) と同じ。
print(A)
print(B)
print("np.tensordot(A,B,axes=(1,0))=")
print(res)
print(resre)
#=============
[OUT]
[[1 2 3]
[4 5 6]]
[[1. 1. 1. 1.]
[1. 1. 1. 1.]
[1. 1. 1. 1.]]
np.tensordot(A,B,axes=(1,0))=
[[ 6. 6. 6. 6.]
[15. 15. 15. 15.]]
[[ 6. 15.]
[ 6. 15.]
[ 6. 15.]
[ 6. 15.]]
＜2つ以上の軸を同時につぶせるテンソル積＞
・AとBの２つの軸で、大きさが同じなら、2軸一緒に（axes=2)つぶせます。Aの後ろ２つの軸と、Bの前２つの軸をつぶすのが、np.tensordot(A, B, axes=2)です。
４次元と３次元で2軸をつぶす（縮約）するテンソル積は4-2+3-2=3次元になります。
３次元と２次元で2軸をつぶす（縮約）するテンソル積は3-2+2-2=1次元になります。
つぶす軸を添え字の番号で指定することもできます。
(例）
# A: 4次元テンソル (2, 3, 4, 5)
# B: 3次元テンソル (4, 5, 2)
A = np.random.rand(2, 3, 4, 5)
B = np.random.rand(4, 5, 2)
# Aの「軸2, 3」(4, 5) と Bの「軸0, 1」(4, 5) を一気にぶつけて消す
res = np.tensordot(A, B, axes=([2, 3], [0, 1])) #axesのあとは添え字の位置の指定
# 結果のテンソルの型は (2, 3, 2) と3次元になります。
print(A)
print(B)
print("np.tensordot(A,B,axes=(1,0))=")
print(res)
[OUT]
略
課題：つぶれる前とつぶれた後のテンソルを視覚化するにはどうしたらよいでしょうか。
たとえば、
aが３軸 (2, 2, 2)のテンソルで成分が、[[[1 2], [3 4]], [[5 6], [7 8]]]だとします。
Aが２軸 (2, 2)のテンソルで成分が、[[1000 100], [10 1]]とするとき、
軸ごとにmatplotlib.pyplotの箱(voxel)を用意すると、
テンソルaは3軸で、2×2×2＝８個の箱（ボクセル）でできた立方体になりますね。
テンソルAは2軸で、2×2＝4個の箱（ボクセル）が1段並びます。
これをツールとして使うことで、テンソル積の軸の選び方で、結果がどう変わるかを確かめましょう。
import numpy as np
import matplotlib.pyplot as plt
import numpy as np
#tensordotの視覚化。axes=1
a = np.array(range(1, 9))
a.shape = (2, 2, 2)
A = np.array((1000,100,10,1), dtype=object)
A.shape = (2, 2)
print("--- a ---")
print(f"形: {a.shape}")
print(a)
print(a[:,-1])
print("\n--- A ---")
print(f"形: {A.shape}")
print(A)
print(A[:,0])
print("\n--- res ---")
fig, ax = plt.subplots(1, 2, figsize=(11, 4), subplot_kw={'projection': '3d'})
def drawTensor3D(P,n):
z_size, y_size, x_size = P.shape[0], P.shape[1],P.shape[2]
voxels = np.ones((x_size, y_size, z_size), dtype=bool)
ax[n].voxels(voxels, edgecolor='k', facecolors='green', alpha=0.1)
#テンソルの内容をその位置にかく
for pz in range(z_size):
for py in range(y_size):
for px in range(x_size):
contents = P[pz,py,px]
ax[n].text(px, py, pz+0.5, f'{contents}', fontsize=15, verticalalignment="bottom")
ax[n].set_xlim([0, 2]);ax[n].set_ylim([0, 2]);ax[n].set_zlim([0, 2])
ax[n].set_xlabel('x');ax[n].set_ylabel('y');ax[n].set_zlabel('z')
ax[n].set_aspect('equal')
def drawTensor2D(P,n):
y_size, x_size = P.shape[0],P.shape[1]
#テンソルの内容をその位置にかく
for py in range(y_size):
for px in range(x_size):
contents = P[py,px]
ax[n].text(px, py, 0.5, f'{contents}', fontsize=15, verticalalignment="top")
ax[n].set_xlim([0, 2]);ax[n].set_ylim([0, 2]);ax[n].set_zlim([0, 2])
ax[n].set_xlabel('x');ax[n].set_ylabel('y');ax[n].set_zlabel('z')
ax[n].set_aspect('equal')
def drawTensor1D(P,n):
x_size = P.shape[0]
#テンソルの内容をその位置にかく
for px in range(x_size):
contents = P[px]
ax[n].text(px, 0, 0, f'{contents}', fontsize=15, verticalalignment="bottom")
ax[n].set_xlim([0, 2]);ax[n].set_ylim([0, 2]);ax[n].set_zlim([0, 2])
ax[n].set_xlabel('x');ax[n].set_ylabel('y');ax[n].set_zlabel('z')
ax[n].set_aspect('equal')
（例１）tensordot(a, A,axes=2)とすると、
aのｚ軸が残り、aのxy軸とAのxy軸がピッタリと重なってつぶれることで、2個の数値になります。
aの[[1 2], [3 4]]がAの[[1000 100], [10 1]]に要素ごとの内積になり
1×1000＋２×100＋3×10＋４×１＝1234ができます。段ｚを１にしてaの2段目から5678ができます。
#============ ここからを変える。================
res = np.tensordot(a, A,axes=2) #aの末２軸とAの前２軸をつぶす。
print(f"形: {res.shape}")
print(res)
drawTensor3D(a,0)
drawTensor2D(A,0)
drawTensor1D(res,1)
plt.show()
#===========================================
--- a ---
形: (2, 2, 2)
[[[1 2]
[3 4]]
[[5 6]
[7 8]]]
[[3 4]
[7 8]]
--- A ---
形: (2, 2)
[[1000 100]
[10 1]]
[1000 10]
--- res ---
形: (2,)
[1234 5678]
（例2）einsum('ijk,lk->ijl',a,A)
einsumはtensordotの添え字番号をijkのように文字にしたもです。
ijkはaの添え字で、lkはAの添え字で、共通の軸はkです。
テンソルaでのkはz軸,y軸,x軸の順のxなので、x軸方法のベクトルとします。
テンソルaは４本のベクトル(1,2),(3,4),(5,6),(7,8)がzｙ平面に立ってるイメージです。
テンソルAでのkはy軸,x軸の順のｘなので、ｘ軸方向の2本のベクトル(1000,100),(10,1)が寝てます。
これらの直積を内積によりつぶすと、１×1000＋２×100＝1200、1×10+2×1＝12、
3400、34、5600、56、7800、78ができますが。これらを2×2×2の箱に入れます。
np.tensordot(a, A, (2, 1))でも同じですね。
#============ ここからを変える。================
res = np.einsum('ijk,lk->ijl',a,A)
print(f"形: {res.shape}")
print(res)
fig, ax = plt.subplots(1, 2, figsize=(11, 4), subplot_kw={'projection': '3d'})
drawTensor3D(a,0)
drawTensor2D(A,0)
drawTensor1D(res,1)
plt.show()
#=========================================
[OUT]
--- a ---
形: (2, 2, 2)
[[[1 2]
[3 4]]
[[5 6]
[7 8]]]
[[3 4]
[7 8]]
--- A ---
形: (2, 2)
[[1000 100]
[10 1]]
[1000 10]
--- res ---
形: (2, 2, 2)
[[[1200 12]
[3400 34]]
[[5600 56]
[7800 78]]]
（例3）einsum('ijk,li->jkl',a,A)
ijkはaの添え字で、liはAの添え字で、共通の軸はiです。
テンソルaでのiはz軸,y軸,x軸の順のzなので、z軸方法のベクトルとします。
テンソルaは４本のベクトル(1,5),(2,6),(3,7),(4,8)がxｙ平面に立っているイメージです。
テンソルAでのiはy軸,x軸の順のｘなので、ｘ軸方向の2本のベクトル(1000,100),(10,1)が寝てます。
これらの直積を内積によりつぶすと、１×1000＋5×100＝1500、1×10+5×1＝15、
3700、37、2600、26、4800、48ができますが。これらを2×2×2の箱に入れます。
np.tensordot(a, A, (0, 1))でも同じですね。
#============ ここからを変える。================
res = np.einsum('ijk,li->jkl',a,A)
print("\n--- res ---")
print(f"形: {res.shape}")
print(res)
drawTensor3D(a,0)
drawTensor2D(A,0)
drawTensor3D(res,1)
plt.show()
#=========================================
[OUT]
--- res ---
形: (2, 2, 2)
[[[1500 15]
[2600 26]]
[[3700 37]
[4800 48]]]
（例４）einsum('ijk,lk->ijl',a,A)
ijkはaの添え字で、liはAの添え字で、共通の軸はkです。
テンソルaでのiはz軸,y軸,x軸の順のzなので、z軸方法のベクトルとします。
テンソルaは４本のベクトル(1,5),(2,6),(3,7),(4,8)がxｙ平面に立っているイメージです。
テンソルAでのiはy軸,x軸の順のｘなので、ｘ軸方向の2本のベクトル(1000,100),(10,1)が寝てます。
これらの直積を内積によりつぶすと、１×1000＋5×100＝1500、1×10+5×1＝15、
3700、37、2600、26、4800、48ができますが。これらを2×2×2の箱に入れます。
np.tensordot(a, A, (0, 1))でも同じですね。
#============ ここからを変える。================
res = np.einsum('ijk,lk->ijl',a,A)
print("\n--- res ---")
print(f"形: {res.shape}")
print(res)
drawTensor3D(a,0)
drawTensor2D(A,0)
drawTensor3D(res,1)
plt.show()
#=========================================
[OUT]
--- res ---
形: (2, 2, 2)
[[[1200 12]
[3400 34]]
[[5600 56]
[7800 78]]]
課題：Geogebraで回転を計算して視覚化するにはどうしますか。
geogebraの数式に、1行ずつx、y、zの関数として入力します。
v1=x*y*z*x;v2=x*y*z*y;v3=x*y*z*z
rotV={Derivative(v3,y)-Derivative(v2,z),Derivative(v1,z)-Derivative(v3,x),Derivative(v2,x)-Derivative(v1,y)}
が回転です。
it=Sequence(t,t,-20,20,10)で、数列を作ります。これを3次元の直積にするために、
st=Flatten(Zip(Zip(Zip((p,q,r),p,it),q,it),r,it))とかきます。p,q,rが独立して動くので直積となり、
テンソルができます。３次元にしたので、stは3軸のテンソルになりますね。
このテンソルは(x,y,z)の集まりともいえますが、
この要素のi番目のｘ座標にアクセスするためにx(st(i))とかきます。
だから、stのi番目の点の座標は(x(st(i)),y(st(i)),z(st(i)))です。この点をPとして、
PからrotVの成分だけ加算した点Qまでつなぐと、ベクトルになりますね。
rotVのｘ成分、ｙ成分、ｚ成分は
Derivative(v3,y)-Derivative(v2,z)、
Derivative(v1,z)-Derivative(v3,x)、
Derivative(v2,x)-Derivative(v1,y)なので、
Qのｘ座標,y座標,z座標は
r1(x,y,z)=x+(Derivative(v3,y)-Derivative(v2,z))/300
r2(x,y,z)=y+(Derivative(v1,z)-Derivative(v3,x))/300
r3(x,y,z)=z+(Derivative(v2,x)-Derivative(v1,y))/300
とおくと、x(st(i)),y(st(i)),z(st(i))をメッシュのテンソルと考えて、ベクトルの配列を作れますね。
それは、
Sequence(Vector((x(st(cnt)),y(st(cnt)),z(st(cnt))),(r1(x(st(cnt)),y(st(cnt)),z(st(cnt))),r2(x(st(cnt)),y(st(cnt)),z(st(cnt))),r3(x(st(cnt)),y(st(cnt)),z(st(cnt))))),cnt,1,a)
です。
とても長ったらしいですが、構造は
Sequence(Vector(P,Q),cnt,1,a)です。
Pをcntの関数として、(x(st(cnt)),y(st(cnt)),z(st(cnt)))とかいて動かし、
Qもcntの関数(r1(x(st(cnt)),y(st(cnt)),z(st(cnt))),r2(x(st(cnt)),y(st(cnt)),z(st(cnt))),r3(x(st(cnt)),y(st(cnt)),z(st(cnt))))
を動かしています。
できなくはないですが、大仕掛けですね。
300で発散の成分を割っているのは大きすぎてはみ出して、動きがわからなくなるからです。
課題：ストークスの定理を視覚化するにはどうしますか。
ガウスの定理、グリーンの定理、ストークスの定理、どれも次元下げの定理です。
だから、空間を扱うときに、空間全部ではなく、表面のようすを調べると、その中は打ち消されているので、全体のようすがわかるというものです。
部分部分をたし込むときに、次元を下げて積分してみましょう。
（例）
from sympy import symbols, cos, sin, Matrix,diff,integrate,pi,lambdify
import numpy as np
import matplotlib.pyplot as plt
#from sympy.plotting import plot3d_parametric_line,plot3d_parametric_surface
# グリーンの定理を３Dにしたものがストークスの定理。曲面のふちだけで積分しても同じ量になる。
# 境界線Cの曲面Sでのベクトル場Fについて、曲面SでのFの回転の面積分G＝Fの各成分の境界線Cで線積分の和
# F=[F1,F2,F3] の回転の面積分。ｎは単位法線ベクトル,Rは曲面Sの領域,Cは領域Rの境界線
# G=∫∫R rotF・ｎ dA = ∫C(F1dx + F2dy + F3dz)=∫C F・dr/ds ds=∫C F・dr
# Sは 放物面 z=f(x,y)=1-(x^2+y^2)(z≧0)で、境界線Cは円周（x^2+y^2=1）F=[y,z,x]
def grad(w,vars): #スカラー場ｗの勾配はベクトル
return Matrix([diff(w, i) for i in vars])
def div(v,vars): #ベクトル場ｖの発散はスカラー
return sum([diff(v[i],vars[i]) for i in range(len(vars))])
def rot(F, vars):#回転はcurlとすることが多い。日本ではrot
F0,F1,F2 = F
return Matrix([diff(F2,vars[1]) - diff(F1,vars[2]),diff(F0,vars[2]) - diff(F2,vars[0]),diff(F1,vars[0]) - diff(F0,vars[1])])
f,xx,yy,zz,x,y,z,r,t,u,v,a,b= symbols('f xx yy zz x y z r t u v a b')
#領域Vで3重積分
F1= y;F2= z;F3= x
F = Matrix([F1,F2,F3])
xx= cos(t)
yy= sin(t)
zz= 0
drdt = Matrix([diff(xx,t),diff(yy,t),diff(zz,t)])
Ft=F.subs({x:xx,y:yy,z:zz})
Ftdotdr=Ft.dot(drdt)
G=integrate(Ftdotdr,(t,0,2*pi))
print("線積分 ∫ F・(dr/dt) dt =",G)
Z=1-(x**2+y**2)
N=grad(z-(Z),[x,y,z])
xxr= r*cos(t)
yyr= r*sin(t)
zzr= 0
rotF = rot(F,[x,y,z])
FN=rotF.dot(N).subs({x:xxr,y:yyr,z:zzr})
FNr=FN*r
G2=integrate(integrate(FNr,(r,0,1)),(t,0,2*pi))
print("面積分 ∫∫ rotF・N dxdy =∫∫ rotF・N(r,t) rdrdt =",G)
#視覚化
def drawVecAt(u0,v0):
#接点と法線の方向
x0, y0, z0 = np.array(func_pos(u0, v0), dtype=float).flatten()
nx, ny, nz = np.array(func_normal(x0, y0), dtype=float).flatten()
fx, fy, fz = np.array(func_F(x0, y0, z0), dtype=float).flatten()
ax.scatter(x0, y0, z0, color= "green", s=10)
#法線ベクトル
ax.quiver(x0, y0, z0, nx/10, ny/10, nz/10, length=1.0, color="green", linewidth=1)
#Fベクトル
ax.quiver(x0, y0, z0, fx, fy, fz, length=1.0, color="blue", linewidth=1)
XXX= r*cos(t)
YYY= r*sin(t)
ZZZ= 1-r**2
N=grad(z-(Z),[x,y,z])
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(111, projection='3d')
t_vals = np.linspace(0, 2*np.pi , 120);
r_vals = np.linspace(0, 1 , 120);
T, R = np.meshgrid(t_vals, r_vals)#２パラメータのテンソル化
func_pos = lambdify((t, r), [XXX,YYY,ZZZ], "numpy")
func_normal = lambdify((x, y), N, "numpy")
func_F = lambdify((x, y, z), F, "numpy")
X, Y, Z = func_pos(T, R)
ax.plot_surface(X, Y, Z, alpha=0.3, cmap='viridis', edgecolor='blue')
pT = np.random.uniform(0, 2*np.pi, size=(150, 1))
pR = np.random.uniform(0, 1.0 , size=(150, 1))
for i in range(len(pT)):
drawVecAt(pT[i].flatten(),pR[i].flatten())
ax.set_xlim([-1.2, 1.2])
ax.set_ylim([-1.2, 1.2])
ax.set_zlim([0, 1])
ax.set_xlabel('x');
ax.set_ylabel('y');
ax.set_zlabel('z')
ax.set_aspect('equal')
ax.set_title("STOKES’S THEOREM")
plt.show()
#=====================================
[OUT]線積分 ∫ F・(dr/dt) dt = -pi
面積分 ∫∫ rotF・N dxdy =∫∫ rotF・N(r,t) rdrdt = -pi