ガウス曲率K、平均曲率Hの計算

Differential Forms with Applications to the Physical Sciencesのp.48の6番でz=f(x,y)と表せるような曲面についてガウス曲率Kと平均曲率Hを求める問題がある。これをsympyでやってみる。

from sympy import *
def cross(X, Y):
    return [X[1]*Y[2]-X[2]*Y[1], X[2]*Y[0]-X[0]*Y[2], X[0]*Y[1]-X[1]*Y[0]]

# 文字の定義
x, y, dx, dy = symbols("x y dx dy")
f = Function("f")(x, y)
fx, fy, fxx, fxy, fyy = f.diff(x), f.diff(y), f.diff(x, x), f.diff(x, y), f.diff(y, y)
# 表示のための置き換えの定義
expsub = {Derivative(f, x):"fx", Derivative(f, y):"fy", Derivative(f, x, x):"fxx", Derivative(f, x, y):"fxy", Derivative(f, y, y):"fyy"}

# v1、v2が曲面の接線、v3が法線
v1 = [i / sqrt(1 + fx**2) for i in [1, 0, fx]]
v2 = [i / sqrt(1 + fy**2) for i in [0, 1, fy]]
v3 = [i / sqrt(1 + fx**2 + fy**2) for i in [-fx, -fy, 1]]
# 局所で正規直行基底を計算
e1, e3 = v1, v3
e2 = cross(v3, v1)
de3 = [i.diff(x) * dx + i.diff(y) * dy for i in e3]

# 接ベクトルを局所座標で表すときの係数を計算
dX = [i*dx + j*dy for i,j in zip([1, 0, fx], [0, 1, fy])]
s1, s2 = Matrix([e1[:2], e2[:2]]).T.inv() * Matrix(dX[:2])
# 計算の確認
if simplify(e1[2]*s1 + e2[2]*s2 - dX[2]) != 0:
   print("wrong calculation")
   import sys
   sys.exit()

# 法線ベクトル
w1, w2 = Matrix([e1[:2], e2[:2]]).T.inv() * Matrix(de3[:2])
if simplify(e1[2]*w1 + e2[2]*w2 - de3[2]) != 0:
   print("wrong calculation")
   sys.exit()

for i in (s1, s2, w1, w2):
    print(simplify(i).subs(expsub))

を実行すれば

(dx*(fx**2 + 1) + dy*fx*fy)/sqrt(fx**2 + 1)
dy*sqrt(fx**2 + fy**2 + 1)/sqrt(fx**2 + 1)
-(dx*fxx + dy*fxy)/(sqrt(fx**2 + 1)*sqrt(fx**2 + fy**2 + 1))
(-dx*fx**2*fxy + dx*fx*fxx*fy - dx*fxy - dy*fx**2*fyy + dy*fx*fxy*fy - dy*fyy)/(sqrt(fx**2 + 1)*(fx**2 + fy**2 + 1))

が得られ、、であることが分かっているので(手計算...)、

$$ \begin{align} K &= &\frac{f_{xx}f_{yy}-f_{xy}^{2}}{(1+f_{x}^{2}+f_{y}^{2})^{2}}\\ H &= &\frac{(1+f_{x}^{2})f_{yy}-2f_{x}f_{y}f_{xy}+(1+f_{y}^{2})f_{xx}}{(1+f_{x}^{2}+f_{y}^{2})^{3/2}} \end{align} $$ が導ける。

参考文献

• SymPy vs. Maxima · sympy/sympy Wiki · GitHub: sympyとmaximaの比較。微分形式の計算とかmaximaで出来ないかな?