2015年1月21日水曜日

wxMaximaで2次元球面のガウス曲率を計算する

前回の計算を発展させて2次元球面のガウス曲率を計算してみる。

本体で読み込むライブラリを作成する

前回と同様、2sphere.macが次の内容だとする。これを「c:\Users\username\Documents\maxima」フォルダ内に置く。

x(\theta,\phi):=r*sin(\theta)*cos(\phi)$
y(\theta,\phi):=r*sin(\theta)*sin(\phi)$
z(\theta,\phi):=r*cos(\theta)$

define(x1(\theta,\phi),diff(x(\theta,\phi),\theta))$
define(y1(\theta,\phi),diff(y(\theta,\phi),\theta))$
define(z1(\theta,\phi),diff(z(\theta,\phi),\theta))$

define(x2(\theta,\phi),diff(x(\theta,\phi),\phi))$
define(y2(\theta,\phi),diff(y(\theta,\phi),\phi))$
define(z2(\theta,\phi),diff(z(\theta,\phi),\phi))$

e1(\theta,\phi):=[x1(\theta,\phi),y1(\theta,\phi),z1(\theta,\phi)]$
e2(\theta,\phi):=[x2(\theta,\phi),y2(\theta,\phi),z2(\theta,\phi)]$

E(\theta,\phi):=trigsimp(e1(\theta,\phi).e1(\theta,\phi))$
F(\theta,\phi):=trigsimp(e1(\theta,\phi).e2(\theta,\phi))$
G(\theta,\phi):=trigsimp(e2(\theta,\phi).e2(\theta,\phi))$



また、3次元ベクトル空間でクロス積とベクトルの大きさを計算できるようにしたい。次の内容のvector3cross.macを作る。これを上のファイルと同じ「c:\Users\username\Documents\maxima」フォルダ内に置く。

cross(a,b):=[a[2]*b[3]-a[3]*b[2],a[3]*b[1]-a[1]*b[3],a[1]*b[2]-a[2]*b[1]]$
norm(v):=sqrt(v.v)$


ガウス曲率を計算する本体を作成する

最後に、ガウス曲率を計算する2sphere2.macを以下の内容で作成する。上で定義した2つのファイルを利用する。

?xchdir("c:/Users/username/Documents/maxima")$
batch("vector3cross.mac")$
batch("2sphere.mac")$

n(\theta,\phi):=trigsimp(cross(e1(\theta,\phi),e2(\theta,\phi)))$
e3(\theta,\phi):=n(\theta,\phi)/(sin(\theta)*r^2)$

define(x11(\theta,\phi),diff(x1(\theta,\phi),\theta))$
define(y11(\theta,\phi),diff(y1(\theta,\phi),\theta))$
define(z11(\theta,\phi),diff(z1(\theta,\phi),\theta))$

define(x12(\theta,\phi),diff(x1(\theta,\phi),\phi))$
define(y12(\theta,\phi),diff(y1(\theta,\phi),\phi))$
define(z12(\theta,\phi),diff(z1(\theta,\phi),\phi))$

define(x22(\theta,\phi),diff(x2(\theta,\phi),\phi))$
define(y22(\theta,\phi),diff(y2(\theta,\phi),\phi))$
define(z22(\theta,\phi),diff(z2(\theta,\phi),\phi))$

e11(\theta,\phi):=[x11(\theta,\phi),y11(\theta,\phi),z11(\theta,\phi)]$
e12(\theta,\phi):=[x12(\theta,\phi),y12(\theta,\phi),z12(\theta,\phi)]$
e22(\theta,\phi):=[x22(\theta,\phi),y22(\theta,\phi),z22(\theta,\phi)]$

L(\theta,\phi):=trigsimp(e11(\theta,\phi).e3(\theta,\phi))$
M(\theta,\phi):=trigsimp(e12(\theta,\phi).e3(\theta,\phi))$
N(\theta,\phi):=trigsimp(e22(\theta,\phi).e3(\theta,\phi))$

K(\theta,\phi):=trigsimp((L(\theta,\phi)*N(\theta,\phi)-M(\theta,\phi)^2)/(E(\theta,\phi)*G(\theta,\phi)-F(\theta,\phi)^2))$


実行と計算結果の確認方法

wxMaximaを起動して、「ファイル(F)」 > 「バッチファイル(B)」から、用意したバッチファイル(2sphere2.mac)を選ぶ。すると一連の計算が実行される。

計算したガウス曲率Kの値を見たいときは、読み込んだ後に
K(\theta,\phi);
と入力すれば(;の後はShift+Enterを入力する)、ガウス曲率の値が表示される。


本体ファイルの説明

?xchdir("c:/Users/username/Documents/maxima")$ は、カレントフォルダを変更する。Lisp で定義されている関数 xchdir を Maxima から呼び出して利用する

batch("vector3cross.mac")$
batch("2sphere.mac")$
は、最初に定義した2つのファイルをバッチファイルとして読み込んで実行する。

0 件のコメント :

コメントを投稿