|
VectorCalculus パッケージ
まず、VectorCalculus パッケージを呼び出します。
|
>
|
restart;
with(VectorCalculus):
interface(showassumed=0);
|
|
|
基本のオブジェクト
|
|
このパッケージで操作する主なオブジェクトは、ベクトルです。座標系や座標名といった、オブジェクトに関する特別な情報は、ベクトル上の属性として保存されています。
|
>
|
v := Vector( [a,b,c] );
|
|
>
|
SetCoordinates( polar );
|
前の例では、各点をそれぞれの座標系内で表現しています。ベクトル場を表現するためには、特別の属性が必要となります。
|
>
|
SetCoordinates( spherical[r,phi,theta] );
|
|
>
|
F := VectorField( <r,0,0> );
|
ベクトル場で通常必要となる座標の名前は、座標系の名前の添え字として保存されていることに注意して下さい。また、出力時の違いにも注意が必要です。ベクトル場の表示を行う場合、e[r] オブジェクトは実際には基底ベクトルで表現され、上線をつけて表示されます。しかし、ベクトルの表示を行う場合には、 e[r] オブジェクトはプレースホルダとなり、上線は使用されません。これらのオブジェクトがデカルト座標系内にのみ存在し、ベクトルとベクトル場の違いが霞んでしまう場合には、この区別が必要です。
上記の出力形式を使用したくない場合には、BasisFormat コマンドを用いて設定の切り替えを行うことが可能です。
BasisFormat コマンドは、いくつかの kernelopts の呼び出しのように、以前の値を返します。そして状態の保存や読み込みがとても簡単に行えます。最後に、MapToBasis コマンドは、種種の座標系間でオブジェクトの変換を可能にしています。
|
>
|
SetCoordinates( polar );
|
|
>
|
MapToBasis( <r,theta>, 'cartesian' );
|
|
>
|
SetCoordinates( spherical );
|
|
>
|
MapToBasis( <r,phi,theta>, 'cartesian' );
|
|
>
|
SetCoordinates( spherical[r,phi,theta] );
|
|
>
|
MapToBasis( VectorField( <r,0,0> ), 'cartesian'[x,y,z] );
|
|
|
|
標準的な操作
|
|
3 つのトップレベルの演算子、`+` (加法), `*` (スカラー乗法) および `.` (内積)、が、このパッケージで作成されるベクトル・属性の組み合わせの操作に重複しています。また、新しい演算子、`&x` (外積 が、このパッケージから出力されます。、可能である場合には、これらの演算子は、標準のベクトル微分演算子 Del ( または Nabla ) として理解するように拡張され、またパッケージから出力されます。
基本算術
|
>
|
SetCoordinates( cartesian );
|
|
>
|
(<r(a+h),s(a+h),t(a+h)> - <r(a),s(a),t(a)>) / h;
|
内積
|
>
|
SetCoordinates( polar );
|
|
>
|
SetCoordinates( cartesian[x,y,z] );
|
|
>
|
Del . VectorField( <x^2,y^2,z^2> );
|
|
>
|
(Del . Del)( f(x,y,z) );
|
|
>
|
L := VectorField( <x,y,z> ) . Del;
|
外積
|
>
|
SetCoordinates( cylindrical );
|
|
>
|
SetCoordinates( cartesian[x,y,z] );
|
|
>
|
Del &x VectorField( <y,-x,z> );
|
|
>
|
L := VectorField( <x,y,z> ) &x Del;
|
Curl, Divergence, Gradient, Laplacian 演算子の微分形式を表示するためには、関数に引数を指定しないで下さい。
|
>
|
SetCoordinates( spherical[r,theta,phi] );
|
|
|
|
曲線
|
|
VectorCalculus パッケージは、曲線に関する標準的な計算を実行するための、多くのルーチンを持っています。
例として、楕円の縮閉線を計算します。曲線の縮閉線は、曲線の曲率の中心に位置するように定義されます。まず、楕円に関していくつかの項目を計算する必要があります。
|
>
|
SetCoordinates( cartesian );
|
|
>
|
ell := <2*cos(t),sin(t)>;
|
|
>
|
nv := simplify( PrincipalNormal(ell,t) );
|
|
>
|
len := simplify( LinearAlgebra:-Norm( nv, 2 ) );
|
|
>
|
r := simplify( RadiusOfCurvature(ell) );
|
ここで、必要な量が得られるので、評価を行うための標準的な公式を使用します。
|
>
|
ev := simplify( ell + r * nv / len );
|
|
>
|
plot( [ [ell[1], ell[2], t=0..2*Pi], [ev[1], ev[2], t=0..2*Pi] ] );
|
R^3 内の曲線が平面内にあるか、捩れなしである場合の、計算の他の応用例を示します。
|
>
|
simplify( Torsion( <cos(t),sin(t),t> ) );
|
|
>
|
simplify( Torsion( <cos(t),sin(t),4-cos(t)-sin(t)> ) );
|
|
|
|
積分
|
|
警告: 明示的に指定されない仮定があると、予期しない結果が生じることがあります。たとえば、球座標で、半径成分は正と仮定されなければなりません。 これは、assume 関数 assume(radius > 0) を用いて行われます。
このパッケージを用いることで、多くの異なる種類の多変数の積分を計算することが可能です。また、積分の標準領域が、特別の構文を用いてサポートされています。
例として、以下のよく知られた結果から始めます。
|
>
|
Int( exp(-x^2), x=-infinity..infinity ) = sqrt(Pi);
|
R^2 全体において exp(-x^2-y^2) を積分すると、Pi が得られることが予測されます。R^3 全体において exp(-x^2-y^2-z^2) を積分すると、Pi^(3/2) が得られることが予測されます。これを、R^3 全体で積分することで、試してみて下さい。
|
>
|
int( exp(-x^2-y^2), [x,y] = Circle( <0,0>, R ) );
|
半径が無限大であると指定することも可能です。
|
>
|
int( exp(-x^2-y^2-z^2), [x,y,z] = Sphere( <0,0,0>, infinity ) );
|
VectorCalculus パッケージを用いて、与えられた曲面を通る、ベクトル場の流量を計算します。
|
>
|
r := 'r'; assume( k::posint );
|
|
>
|
F := VectorField( <r^(-k),0,0>, spherical[r,phi,theta] );
|
|
>
|
Flux( F, Sphere( <0,0,0>, R ) ) assuming R>0;
|
別の例として、2 つの立体を場の中心で分離するために必要な作業を計算します。
|
>
|
r := 'r'; k := 'k'; assume(k>0);
|
|
>
|
F := VectorField( <-k/r^2,0,0>, spherical[r,phi,theta] );
|
評価前の不活性形の積分は、関数 ArcLength, Flux, int, LineInt, PathInt, SurfaceInt に対して、パラメータ 'inert' を指定することにより表示できます。
|
>
|
SetCoordinates('cartesian'[x,y,z]);
|
|
>
|
int( exp(-x^2)*ln(x), x=0..infinity, 'inert');
|
|
>
|
LineInt( VectorField( <y,-x,z> ), Circle3D( <0,0,0>, r, <1,1,1> ), 'inert' );
|
|
>
|
F := MapToBasis( F, cartesian[x,y,z] );
|
簡単化のため、積分経路を正の x 軸に沿った直線にとります。
|
>
|
LineInt( F, Line( <eps,0,0>, <r,0,0> ) ) assuming eps>0;
|
評価前の不活性形の積分は、関数 ArcLength, Flux, int, LineInt, PathInt, SurfaceIntに対して、パラメータ 'inert'を指定することにより表示できます。
|
>
|
SetCoordinates('cartesian'[x,y,z]);
|
|
>
|
int( exp(-x^2)*ln(x), x=0..infinity, 'inert');
|
|
>
|
LineInt( VectorField( <y,-x,z> ), Circle3D( <0,0,0>, r, <1,1,1> ), 'inert' );
|
|
|
|
最大化および最小化
|
|
このパッケージのルーチンと LinearAlgebra パッケージは、ラグランジュ乗数の手法を用いた場合の、最大化および最小化の問題を学習するのに役立ちます。例として、R^3 内のある点とある曲面の間の最短距離に関する公式を導き出します。この問題を、拘束つきの最小化問題として設定することが可能です。
|
>
|
SetCoordinates( cartesian[x,y,z] );
|
この問題における点が (p1, p2, p3) で与えられるものと仮定します。すると、最小化したい関数は以下のようになります。
|
>
|
f := (x-p1)^2 + (y-p2)^2 + (z-p3)^2;
|
次に、(x, y, z) が与えられた平面上にあるという条件は次の拘束を与えます:
|
>
|
cond := a*x + b*y + c*z + d;
|
ラグランジュ乗数の手法は、これらの勾配が比例するという指定を行います。
|
>
|
gcond := Gradient( cond );
|
そこで、解くべき連立方程式は以下のようになります:
|
>
|
eqns := { seq( gf[i] = lambda * gcond[i], i=1..3 ), cond = 0 };
|
|
>
|
sols := convert( solve( eqns, {x,y,z,lambda} ), list );
|
ここで、a^2+b^2+c^2 はゼロではないことに注意して下さい。この場合、曲面が存在しないあるいは退化した曲面のいずれかであるためです。次に、この解が最小、最大、またはそのどちらでもないかを判定します。この判定には、ヘッセの判定法を使用します。
|
>
|
h := Hessian( f - lambda*cond, [x,y,z,lambda] );
|
|
>
|
LinearAlgebra:-Determinant( eval( h, sols ) );
|
この行列式は常に負であるため、最小値が求まりました。そこで、この解を元の関数 f に代入して、公式を導きます。
|
>
|
simplify( sqrt( eval( f, sols ) ) ) assuming real;
|
|
|
|
拡張性
|
|
組み込まれている座標系へアクセスするのと同じように、パッケージに座標系を追加することが可能です。これにより、全ての VectorCalculus コマンドが、追加の座標系で計算を行うことが可能となります。必要となるのは、この新しい座標系の単位基底ベクトルが他のものと直交している必要があるという条件だけです。
|
>
|
u := 'u': v := 'v': f := 'f': assume( u>0, v>0 );
|
|
>
|
AddCoordinates( mycoords[u,v], [u^2+v^2,u^2-v^2] );
|
|
>
|
SetCoordinates( mycoords[u,v] );
|
|
>
|
MapToBasis( <u,v>, 'cartesian' );
|
|
>
|
SetCoordinates( cartesian );
|
|
>
|
MapToBasis( <1,1>, 'mycoords' );
|
|
|
