行列計算における数式処理ソフト maxima の利用について
足立健朗

Contents
1. イントロダクション
1.1. maxima とは何か
1.2. maxima の入手とインストール
1.3. maxima の実行と処理、そして終了
1.4. maxima についての詳しい情報
2. 行列を与える
3. 行列の和と積
3.1. 和、スカラー倍
3.2. 行列の積3.3. 行列の巾乗
4. 逆行列の計算
5. 連立一次方程式
5.1. 連立一次方程式を解く
5.2. 行列表現
5.3. 行列の階数と連立一次方程式の解
5.4. 同次連立一次方程式
6. 正方行列6.1. 行列式とトレース
6.2. 転置行列、対称行列、交代行列、エルミート行列、歪エルミート行列
7. maxima 内部における行列の扱いについて
7.1. 配列とリストと行列の関係について
7.2. 行と列の操作について
7.3. 成分の操作について
8. 線形写像と基底の変換
8.1. 線形写像、線形変換
8.2. 核、像、次元
8.3. 基底と成分の変換
9. ベクトルの内積
9.1. ベクトルの内積について
9.2. Gram-Schmidt の直交化
9.3. 直交行列とユニタリ行列
10. 固有値と固有ベクトル10.1. 固有値と固有方程式
10.2. 対称行列の対角化
10.3. 行列の巾乗再論
1

2
2
2
3
5
7
10
10
11
12
14
16
16
16
17
20
22
22
26
30
30
34
42
45
45
48
53
59
59
68
74
82
82
86
93

2

足立健朗

11. ベクトル解析
11.1. 外積
11.2. 微分作用素
12. あとがきというか言い訳

95
95
99
105

行列計算における数式処理ソフト maxima の利用について

3

1. イントロダクション
1.1. maxima とは何か. maxima は汎用の数式処理ソフトである。そして、数式処理ソ
フトとは、数式のまま計算をすることのできる計算機プログラムのことである。電卓的に数式を入力計算できるだけでなく、集合、ベクトル、配列、数列、級数、行列、関数などの様々なデータ型を扱うことができ、それらをプログラムして処理することのできる超高水準プログラミング言語としての機能をもったソフトである。
maxima は MIT(マサチューセッツ工科大学) の Project MAC で開発された Macsymaの子孫である。Macsyma とは MAC’s SYmbolic MAnipulator の略であり、記号積分をコンピュータで実行させる研究の成果として、同じく MIT で開発された人工知能用プロ
グラミング言語 Lisp の応用ソフトの一つとして生まれた数式処理ソフトである。Project
MAC の解散とともに、Macsyma はいくつかの枝に分岐する。その一つは Symbolics Inc.に売却された Macsyma であり、もう一つは米国エネルギー省の管理の下に置かれた DOE
Macsyma である。Symbolics の Macsyma は Symbolics Inc. とその後継者の Macsyma
Inc が 1999 年になくなるまで、商用の数式処理ソフトして販売されていた。Texas 大学の William Schelter は個人的に Doe Macsyma のソースコードを所持し、それを GNU
Common Lisp に移植しメンテナンスしていたが、1990 年代半ばに、それを maxima と名付け公開する決意をした。そして MIT や国防省とのライセンスについての合意の下にそれを公開した。Bill Schelter は 2001 年に、心臓発作のためなくなったが、現在は UCBの R.Fateman らを中心とする、インターネット上のボランティアたちの手によってオー
プンソースの GPL ソフトウェアとして維持、開発されている。
maxima と同様の機能をもった数式処理ソフトは他にもいくつも存在する。主なものをあげると、Asir、Maple、Mathematica、MuPAD、Reduce、Yacas などがある。このうち、Maple、Mathematica は有名な商用ソフトであり、Reduce は Macsyma と同じ時期に開発された商用のソフトである。Asir は日本で開発された数式処理ソフトであるが、得意な分野が限定されているので、他のものほどには汎用性は無い。Yacas は最近開発された GPL ライセンスで配布されるフリーな数式処理ソフトであるが、maxima
ほどの機能はもっていないようである。さて、この文書で特に maxima を採用しているのは、フリーでオープンソースな汎用数式処理ソフトで、ライセンスを気にせず、高額な投資をすることなく教育目的に利用できるという理由からである。したがって、Asir、
MuPAD、Yacas なども同じ目的に使用することができる。
1.2. maxima の入手とインストール. maxima は Linux、 FreeBSD、MacOS X など、
Unix と同等の機能をもつ OS に Common Lisp 処理形がインストールされた環境が存在するならば、ソースコードをダウンロードして、コンパイル、インストールすることができる。もし、あなたが Windows ファミリーの OS しか使えないならば、maxima の実行形式版をダウンロードしてインストールすることによって、maxima を使うことができる。あるいは、Windows に cygwin という、疑似 Unix 環境を提供するソフトをインストールしておいて Unix 版の maxima をインストールすることも可能である。Windows版 maxima 及びソースコード版 maxima は、どちらもhttp://maxima.sourceforge.net/
からダウンロードできる。好みのウェッブブラウザーを使って上記 URL にアクセスしてみるとよい。Windows 版の場合、ダウンロードしたファイルをダブルクリックすればイ
ンストーラが起動するので、適宜質問に答えながらボタンをクリックしていけばインストールできるはずである。インストールが終了すれば、スタートメニューに maxima が

4

足立健朗

登録されているはずである、また、デスクトップ上に maxima のアイコンが作られているはずである。

1.3. maxima の実行と処理、そして終了. Unix 版で Terminal のコマンドラインから
maxima を立ち上げると、次のように表示される。
i i i i i i i
I I I I I I I
I \ ‘+’ / I\ ‘-+-’ /‘-__|__-’|------+------

ooooo
8
8
8
8
8
8
o
ooooo

o
8
8
8
8
8
8oooooo

ooooooo88888
ooo8ooo

ooooo
8
o
8
ooooo
8
o
8
ooooo

ooooo
8
8
8
8
8oooo
8
8
8

Copyright (c) Bruno Haible, Michael Stoll 1992, 1993
Copyright (c) Bruno Haible, Marcus Daniels 1994-1997
Copyright (c) Bruno Haible, Pierpaolo Bernardi, Sam Steingold 1998
Copyright (c) Bruno Haible, Sam Steingold 1999-2002
Maxima 5.9.0 http://maxima.sourceforge.net
Distributed under the GNU Public License. See the file COPYING.
Dedicated to the memory of William Schelter.
This is a development version of Maxima. The function bug_report()
provides bug reporting information.
(C1)
また、Unix 版の xmaxima あるいは Windows 版の maxima を起動すると、専用のウィンドウが開いて、単に(C1)と表示される。いずれの場合も (C1) というのがプロンプト (入力促進記号) であり、これにつづけてあなたのコマンドをタイプ入力していくことにより、対話的に maxima を利用できるようになる。Unix 版の maxima を、TeXmacs や、Emacs+imaxima などの環境で使用すると美しい数式表示が可能になる。以下では、TeXmacs 上で maxima を使用した場合の例をあげることにする。青字の部分がユーザーの入力部分である。
(C1) integrate(1/(1+x^3),x);şť√arctan 2x−1log(x2 −x+1)
3√(D1) −++ log(x+1)
6
3
3(C2) a : (1 + x - 3*y)^5;(D2) (−3y + x + 1)5(C3) expand(D2);(D3) − 243y 5 + 405xy 4 + 405y 4 − 270x2 y 3 − 540xy 3 − 270y 3 + 90x3 y 2 + 270x2 y 2 +270xy 2 + 90y 2 − 15x4 y − 60x3 y − 90x2 y − 60xy − 15y + x5 + 5x4 + 10x3 + 10x2 + 5x + 1(C4) 10!;(D4) 3628800(C5) factor(D4);(D5) 28 34 52 7行列計算における数式処理ソフト maxima の利用について

5

ここで、(C*) につづく文字列が入力したコマンドであり、(D*) に続く文字列が maximaの出力をあらわしている。さらに、上の例では、maxima が履歴を憶えているシステム
であることを利用している。例えば、(C3) における expand(d2); の d2 は (D2) の結果を意味でする。すべての作業が終えるときには、(C*) quit();とする。
maxima を使っていると、しばしばトラブルに出会うだろう。何が起こるかは、状況による。一番目にする可能性が高いのは、あなたがタイプミスをしたときに、表示される
エラーメッセージだ。

(C2) B:matrix([1,2:});
Incorrect syntax: Missing ]B:matrix([1,2:})^(C2) Incorrect syntax: Premature termination of input at ;.
;^こういうときには、慌てずもういちど正しくタイプしなおせばいい。。

(C2) B:matrix([1,2]);
(D2)

[ 1

2 ]

タイプミスでなく、論理的なミスをすることもある。異なったデータ型の変数を足そうとして、(C3) A+B;
FULLMAP found arguments with incompatible structure.
-- an error. Quitting. To debug this try DEBUGMODE(TRUE);)というようなメッセージがあらわれたら、それは入力した式に論理的な誤りがある可能性を示している。入力した式の意味をよく考え直してみよう。もっと深刻なエラーが起きた
場合には、次のようなメッセージがあらわれるかも知れない。

*** - POSITION: :START = 1 should not be greater than :END = 0
The following restarts are available:
R1 = Macsyma top-level
1. Break [1]>これは、maxima では解決できないようなエラーが生じた場合に表示される。

1. Break [1]> R1
あるいは

1. Break [1]> (run)とすることによって maxima がリスタートする。そのまま、maxima を終るには、

1. Break [1]> (bye)
としよう。

6

足立健朗

1.4. maxima についての詳しい情報. maxima を使っていて、しばしば特定のコマンドについてより詳しい情報が欲しくなることがあるだろう。たとえば、コマンド (maxima の正式な用語では関数と言うが、数学の関数とまぎらわしいので、本文書ではコマンドと言うことにする) factor について知りたいとおもったら、プロンプトに続いて describe(“factor”)とタイプする。すると(C1) describe("factor");
0: DONTFACTOR :(maxima.info)Definitions for Matrices and Linear Algebra.
1: EXPANDWRT_FACTORED :Definitions for Simplification.
2: FACTOR :Definitions for Polynomials.
3: FACTORFACSUM :Definitions for Simplification.
4: FACTORFLAG :Definitions for Polynomials.
5: FACTORIAL :Definitions for Number Theory.
6: FACTOROUT :Definitions for Polynomials.
7: FACTORSUM :Definitions for Polynomials.
8: GCFACTOR :Definitions for Polynomials.
9: GFACTOR :Definitions for Polynomials.
10: GFACTORSUM :Definitions for Polynomials.
11: MINFACTORIAL :Definitions for Number Theory.
12: NEXTLAYERFACTOR :Definitions for Simplification.
13: NUMFACTOR :Definitions for Special Functions.
14: SAVEFACTORS :Definitions for Polynomials.
15: SCALEFACTORS :Definitions for Miscellaneous Options.
16: SOLVEFACTORS :Definitions for Equations.
Enter n, all, none, or multiple choices eg 1 3 :のように表示される。(残念ながら英語である。) さらに、自分に必要な情報の番号を選んで、Enter n, all, none, or multiple choices eg 1 3 : 5と入力すれば、
Info from file /usr/local/info/maxima.info:- Function: FACTORIAL (X)The factorial function. FACTORIAL(X) = X! . See also
MINFACTORIAL and FACTCOMB. The factorial operator is !, and the
double factorial operator is !!.
(D1)
FALSE
のように、factor の説明と使い方、場合によっては例などが表示される。
より詳しい情報が知りたいときは、maxima の配布物に含まれる文書、tutorial、manual、info ファイルなどを参照することになる。xmaxima や Windows 版 maxima などを使っていれば, これらをオンラインで参照できる。PDF でフォーマットされた文書が欲しい場合には、以下の URL からダウンロードすればいい。
http://maxima.sourceforge.net/docs.shtml
ここで、Maxima Reference Manual、DOE-Maxima Refernce、The Maxima Book などが手に入る。また、

行列計算における数式処理ソフト maxima の利用について

http://www.bekkoame.ne.jp/~ponpoko/index.htmlには Maxima Reference Manual (maxima-5.6) の日本語訳が置いてある。

7

8

足立健朗

2. 行列を与える
ここでは、maxima において行列を生成する方法について述べる。第 1 の方法は、行列の全ての成分をあからさまに与える方法である。一般的に、行列
a11 a12 · · · a1n
 a21 a22 · · · a2n 
A=
..
.. 
...
 ...
.
. 

am1 am2 · · ·amn

を maxima に与えるには、コマンドラインから A : matrix([a11 , · · · , a1n ],[a21 , · · · , a2n ],· · · ,[am1 , · · · , amn ]); とタイプする。
例 1. (C1)µ A : ¶matrix([a, b], [c, d]);
a b
(D1)
c d
例 2. (C2)2, 3], [0, 4, -1]);µ B : matrix([-1,¶−1 2 3(D2)
0 4 −1

maxima において、コロン ‘“:” は、代入あるいは定義を意味する記号である、コロンの右辺の式にコロンの左辺にある名前をつける、あるいは右辺を左辺の変数の値とするときに使う。
第 2 の方法は行列要素を 2 次元配列と見て、配列から行列を生成する方法である。すな
わち、aij = a[i, j] を与えられた 2 次元配列とするとき、maxima のコマンドラインから
A : genmatrix(a, i2, j2, i1, j1); とタイプすれば行列
ai1,j1 · · · ai1,j2.. ..
A =  ...
.
.
ai2,j1 · · ·ai2,j2が生成される。特に、i1=j1= 1 の場合には、簡略化した構文 A:genmatrix(a, i2, j2); が使える。
例 3. (C1) a[i,j] := i + j - 1;
(D1) ai,j := i + j − 1(C2) 
A : genmatrix(a,3, 4);
1 2 3 4
(D2)  2 3 4 5 
3 4 5 6
注意 1. 配列定義には構文 a[i,j] := “i,j の式” を使う。
注意 2. 配列については 7 節を参照せよ。
第 3 の方法は、maxima に用意されている対話的に行列成分を入力するための関数
entermatrix を用いる方法である。実際例で説明しよう。
例 4. (C1) A : entermatrix(2,3);
Row 1 Column 1: 1;
Row 1 Column 2: 0;

行列計算における数式処理ソフト maxima の利用について

9

Row 1 Column 3: -1;
Row 2 Column 1: 4;
Row 2 Column 2: 2;
Row 2 Column 3: 8;
Matrixµ entered. ¶1 0 −1(D1)
4 2 8
特に正方行列のときには、入力を簡易なものにするための、オプションを選択できる。
例 5. (C2) B : entermatrix(3,3);
Is the matrix 1. Diagonal 2. Symmetric 3.
Answer 1, 2, 3 or 4 : 1;
Row 1 Column 1: -1;
Row 2 Column 2: 2;
Row 3 Column 3: 4;
Matrixentered. −1 0 0(D2)  0 2 0 
0 0 4

Antisymmetric 4.

General

ここでは、B を対角行列として、対角成分だけを入力している。
すでに存在している行列からそのコピーをつくるには、単に代入文を使えばいい。例
えば、(C1) µA : matrix([1,2], [3, 4]);
¶
1 2
(D1)
3 4
(C2) µB : A; ¶
1 2
(D2)
3 4
(C3) µ
B;
¶
1 2
(D3)
3 4
特別な行列を定義するためには、専用のコマンドが用意されている。対角成分がすべて等しい対角行列を定義するには diagmatrix を使う。
例 6. (C1)2);
 A : diagmatrix(4,
2 0 0 0
 0 2 0 0 (D1) 
 0 0 2 0 
0 0 0 2
特に n 次単位行列には ident を用いて、ident(n); とタイプする。
例 7. (C1) ident(3);
1 0 0(D1)  0 1 0 
0 0 1

10

足立健朗

(m, n)-型の零行列を作りたいなら、zeromatrix というコマンドを用いて、A : zeromatrix(m, n); のようにタイプする。
例 8. (C1)zeromatrix(3, 2);
 C :
0 0(D1)  0 0 
0 0
第 (i, j)-成分が与えられた値 a で、他の成分がすべてゼロであるような (m, n) 型行列を与えるコマンドは ematrix(m, n, a, i, j) である。
例 9. (C1)
4, -1, 2, 3);
 D : ematrix(4,
0 0 0 0
 0 0 −1 0 (D1) 
 0 0 0 0 
0 0 0 0(m, 1) 型の行列を m 次元列ベクトルとも言うが、m 次元列ベクトルを与えるコマンドはcolumnvector である。columnvector を使うためには、事前に load(eigen); とタイプして必要なパッケージを読み込ませておく必要がある。
例 10. (C1) load(eigen);
Warning - you are redefining the MACSYMA function EIGENVALUES
Warning - you are redefining the MACSYMA function EIGENVECTORS(D1) /usr/local/share/maxima/5.9.0/share/matrix/eigen.mac(C2) 
a : columnvector([1, 2, 3]);1(D2)  2 
3(C3) x : [1, 2, 3, 4];(D3) [1, 2, 3, 4](C4) 
b : columnvector(x);
1
 2 

(D4) 
 3 
4

行列計算における数式処理ソフト maxima の利用について

11

3. 行列の和と積
3.1. 和、スカラー倍. 同じ型の２つの行列 A, B に対してその和 A + B が定義できる。結果の行列 A + B は A, B と同じ型の行列である。すなわち
a11 a12
 a21 a22

·
 ·
A=
·
 ·
 ·
·
am1 am2
のとき、

···
···
aij
···


a1n
a2n 

· 
,
· 
· 
amn



b11 b12 · · · b21 b22 · · ·

·
 ·
B=
·
bij
 ·
 ·
·
bm1 bm2 · · ·


b1n
b2n 

· 

· 
· 
bmn


a11 + b11 a12 + b12···
a1n + b1n···
a2n + b2n  a21 + b21 a22 + b22


·
·
·


A+B =

·
·
a
+
b
·


ij
ij


·
·
·
am1 + bm1 am2 + bm2···
amn + bmn
である。
maxima で行列の計算をするのは非常に簡単である。A と B をすでに定義された行列とするとき、行列 C = A + B を求めるには、プロンプトから C : A + B; とタイプすればいい。
例 11. (C1)µ A :¶ matrix([a, b], [c, d]);
a b
(D1)
c d
(C2) µB : matrix([1,2], [3, 4]);
¶
1 2
(D2)
3 4
(C3) µC : A + B;¶
a+1 b+2(D3)
c+3 d+4
行列のスカラー倍も簡単である、A を任意の行列、c をスカラーとするとき cA を求めるには、c*A; とタイプするだけでいい。上の例から引続き、例 12. (C4)µ 3*A; ¶
3a 3b
(D4)
3c 3d
(C5) µ
f*A;
¶
af bf
(D5)
cf df
(C6) µ
n*B;
¶
n 2n
(D6)
3n 4n

12

足立健朗

3.2. 行列の積. ２つの行列の積はちょっと難しい。勝手な 2 つの行列に対しては積が定義できるとは限らない。積の定義はつぎの通りである。A = (ai ) を (m, l) 型の行列、B = (bij )を (l, n) 型の行列とするとき、A と B の積 AB = (cij ) が定義できて、AB は (m, n) 型の行列になる。各成分 cij は公式
(1)

cij =

l
X

aik bkj

k=1(2)= ai1 b1j + ai2 b2j + · · · + ail blj

で与えられる。注意すべきことは、AB が定義されても BA が定義されるとは限らない
し、また定義される場合でも AB = BA がなりたつとは限らないことである。
特別な場合として A,B がともに n 次正方行列である場合には、AB,BA はともに定義可能で、ともに n 次正方行列になるが、この場合でも、AB = BA とは限らない。このよ
うな事情を行列の積は非可換であるという。
行列の積を maxima に計算させるには “.” 作用素を使う。A と B が与えられていてAB が計算可能であるとき、maxima では、A.B とタイプする。AB が計算可能でないときはエラーとなる。
例をあげよう。
例 13. (C1)3, 4], [-1, 0, 2]);
µ A : matrix([2,
¶
2 3 4
(D1)
−1 0 2(C2) 
B : matrix([3,0], [-1, 2], [3, 4]);
3 0(D2)  −1 2 
3 4(C3) µC : A.B; ¶
15 22
(D3)
3 8
(C4) 
D : B.A;69 12(D4)  −4 −3 0 
2
9 20(C5) C.D;
incompatible dimensions - cannot multiply-- an error. Quitting. To debug this try DEBUGMODE(TRUE);)(C6) µE : matrix([a,b], [c, d]);
¶
a b
(D6)
c d
(C7) µF : matrix([e,f], [g, h]);
¶
e f
(D7)
g h
(C8) µE.F;¶
bg + ae bh + af
(D8)
dg + ce dh + cf
(C9) F.E;

行列計算における数式処理ソフト maxima の利用についてµ(D9)

cf + ae df + be
ch + ag dh + bg

13¶3.3. 行列の巾乗. 行列の積の定義から、もし行列 A が n 次正方行列であれば、A の k 乗を、帰納的に Ak := AAk−1 で定義できる。maxima はこのような行列の巾乗を計算する手段も備えている。
与えられた正方行列 A に対して、Ak を求めるには Aˆˆk とタイプすればいい。ただし、k は具体的な自然数でなくてはならない。また、Aˆ k とすると各成分を k 乗した行列が得られるがこれは Ak ではないことに注意する。
例 14. (C1)matrix([2, 3], [0, -1]);µ A : ¶
2 3
(D1)
0 −1
(C2) µ
A^^5;
¶
32 33
(D2)
0 −1
(C3) 
B : matrix([a,b, c], [d, e, f], [g, h, i]);a b c(D3)  d e f 
g h i(C4) 
B^3;

a3 b3 c3
(D4)  d3 e3 f 3 
g 3 h3 i3
B 3 は紙に収まりきらないので、端末版の表示をあげよう。

(C6) B^^3;
(D6)
[
2 ]
[ g (c i + b f + a c) + d (c h + b e + a b) + a (c g + b d + a ) ][]Col 1 = [2]
[ g (f i + e f + c d) + d (f h + e + b d) + a (f g + d e + a d) ]
[
]
[
2
]
[ g (i + f h + c g) + d (h i + e h + b g) + a (g i + d h + a g) ][2 ][ h (c i + b f + a c) + e (c h + b e + a b) + b (c g + b d + a ) ][]Col 2 = [2]
[ h (f i + e f + c d) + e (f h + e + b d) + b (f g + d e + a d) ]
[
]
[
2
]
[ h (i + f h + c g) + e (h i + e h + b g) + b (g i + d h + a g) ]

14

足立健朗

[
2 ]
[ i (c i + b f + a c) + f (c h + b e + a b) + c (c g + b d + a ) ][]Col 3 = [2]
[ i (f i + e f + c d) + f (f h + e + b d) + c (f g + d e + a d) ]
[
]
[
2
]
[ i (i + f h + c g) + f (h i + e h + b g) + c (g i + d h + a g) ]これらの計算は機械的に行なわれるので、大きな行列の巾乗とか、巾乗の指数が大きい
場合には計算機のリソースを大量に消費する。したがって、安易に使うことは推奨できない。理論的なアプローチ、例えば座標変換して対角行列にするなどのアプローチが可能かどうかを調べるなどすべきである。

行列計算における数式処理ソフト maxima の利用について

15

4. 逆行列の計算
正方行列 A に対して、AX = XA = I を満たす A と同じ型の正方行列 X のことを、
A の逆行列といって A−1 と書く、ただし I は単位行列である (高校では E と書くかもしれない)。逆行列を maxima で形式的に求めることは簡単にできる。A を与えられた正方行列とするとき、A−1 を maxima で求めるには、Aˆˆ(-1); あるいは invert(A); とタイプすればよい。
例 15. (C1)µ A :¶ matrix([a, b], [c, d]);
a b
(D1)
c d
(C2) µA^^(-1);¶db− ad−bc
ad−bc
(D2)
c
a
− ad−bcad−bc(C3) µB : matrix([1,2], [2, 1]);
¶
1 2
(D3)
2 1
(C4) µB^^(-1); ¶− 13 23
(D4)
2
− 13
3
(C5) µinvert(B);¶− 13 23
(D5)
2
− 13
3
ここで注意したいのは、もし ∆ = ad − bc = 0 ならば A−1 は存在しないのであるが、maxima は気にせず形式的計算を実行していることである。このあたりが、真の数学と計算機の理解している数学の差であり、使うものはよく注意する必要がある。これらは良く知られた Cramer の公式を用いて計算されている。これまた良く知られているように、 Cramerの公式は理論的には重要であるが、計算量が大きいので、巾乗のところで注意したことがこの場合にもあてはまる。すなわち、大きな行列に対して安易にこの方法で計算することは勧められない。行列の基本変形などを使うことを推奨する。さて、付録として Cramerの公式を maxima で表現しておこう。expand(adjoint(A))/expand(determinant(A)) により A の逆行列を計算するのが Cramer の公式である。ここで adjoint(A) は A の余因子
行列を、determinant(A) は A の行列式を表す。expand は多項式を展開するコマンドである。
例 16. (C1)µ A :¶ matrix([a, b], [c, d]);
a b
(D1)
c d
(C2) µexpand(adjoint(A))/expand(determinant(A));¶
db− ad−bc
ad−bc
(D2)
a
c
− ad−bcad−bc(C3) 
C : matrix([a,b, c], [d, e, f], [g, h, i]);a b c(D3)  d e f 
g h i(C4) ratsimp( invert(C) - expand(adjoint(C))/expand(determinant(C)) );

16

足立健朗




0 0 0
(D6)  0 0 0 
0 0 0
ここで、最後に使ったコマンド ratsimp は有理多項式を計算して簡約表示するものである。余因子行列、行列式については後で再び取り上げる。

行列計算における数式処理ソフト maxima の利用について

17

5. 連立一次方程式
5.1. 連立一次方程式を解く. n 個の未知数 x1 , x2 , . . . , xn について m 個の一次方程式からなる連立一次方程式
a11 x1 + a12 x2 + · · · + a1n xn
= b1



a x + a x + · · · + a x
= b2
21 1
22 2
2n n
(3)

······
am1 x1 + am2 x2 + · · · + amn xn = bm
を maxima で取り扱う方法について説明しよう。(3) を解くにはコマンド solve を使う。
solve は汎用の方程式を解くためのコマンドで、連立ー次方程式に限らず使用することができる。(解が求まるとは限らないが！) 一般的な構文は(C1) solve([方程式１, 方程式 2, ... , 方程式 m], [未知数１, 未知数２, ... , 未
知数 n]);である。ここで、方程式はそのままタイプすると面倒なので、eq1 : 方程式１ のように名前をつけて、solve([eq1, eq2, ... , eq2], [x, y, ... , z]); のように使うといい。
例 17. (C1) eq1 : a*x + b*y = 0;
(D1) by + ax = 0(C2) eq2 : c*x + d*y = 0;
(D2) dy + cx = 0(C3) solve([eq1, eq2], [x, y]);(D3) [[x = 0, y = 0]](C4) eq3 : a*x + b*y = e;(D4) by + ax = e(C5) eq4 : c*x + d*y = f;(D5) dy + cx = f(C6) ££
solve([eq3, eq4], [x,¤¤ y]);
de−bf
ce−af
(D6)
x = − bc−ad , y = bc−ad(C7) ££
solve([3*x-2*y+z=3,¤¤x+y-3*z=7, y+z=9], [x, y, z]);(D7)
x = 26, y = 36
, z = 95
55(C8) ££
solve([3*x-2*y+z=3,¤¤x+y-3*z=7], [x, y]);(D8)
x = 5z+17, y = 10z+18
55(C9) solve([3*x - 2*y = 3, x + y = 7, x - y = 9], [x, y]);Inconsistent equations: (3)-- an error. Quitting. To debug this try DEBUGMODE(TRUE);)
上の (C8) は未知数の数が方程式の数より多い場合であるが、このような場合には一般的に解に任意パラメータが含まれる。ここでは、z が任意パラメータである。逆に (C9) は
未知数の数よりも方程式の数の方が多い過剰決定系とよばれる場合であるが、多くの場合解を持たない。この場合も解をもたないので、maxima が警告を発している。

5.2. 行列表現. 線形代数学で連立一次方程式を理論的に扱う場合には、連立一次方程式を行列とベクトルの間の方程式に表現して扱うのがスマートである。連立一次方程式 (3) に

18

足立健朗

対して、



a11
 a21
A=
 ...

a12
a22
..
.

···
···
...

am1 am2 · · ·


a1n
a2n 
,
.. 
. 




x1
 x2 
x=
 ...  ,

amn

xn




b1
 b2 
b=
 ... bmと置くと、(3) は(4)

Ax = bと書くことができる。これを (3) の行列表現と言う。ここで、A を (3) の係数行列という。
また、次の行列
a11 a12 · · · a1n b1
 a21 a22 · · · a2n b2 (A|b) := 
..
..
.. 
..
 ...
.
.
.
. 

am1 am2 · · ·amn bmを (3) の拡大係数行列という。maxima には、与えられた連立一次方程式から、係数行列および拡大係数行列を取り出すコマンドが用意されている。未知数 x, y, ... , z に対する方程式 eq1, eq2, ... , eqm からなる連立一次方程式の係数行列を求めるための構文は(C1) A : coefmatrix([eq1, eq2, ... , eqm], [x, y, ... , z]);である。これで記号 A は係数行列を意味することになる。また、拡大係数行列を求める
構文は(C1) Ab : augcoefmatrix([eq1, eq2, ... , eqm], [x, y, ... , z]);である。これで記号 Ab は拡大係数行列を意味することになる。
注意 3. maxima の augcoefmatrix は (A|b) ではなく、(A| − b) を求める。しかし、階数の問題については同等である。

5.3. 行列の階数と連立一次方程式の解. 与えられた連立一次方程式が解けるかどうかを、
solve コマンドで実験する以前に調べるために、準備として行列の階数について説明しよう。任意の行列
a11 a12 · · · a1n
 a21 a22 · · · a2n 
A=
,
..
.. 
...
 ...
.
. 
am1 am2 · · ·amnに対して、(1) ある行と別のある行をいれかえる。(2) ある行を定数倍する。(3) ある行に別のある行を加える。という 3 種類の操作を施すことを行列 A の行基本変形という。これらの操作を繰り返すことにより、A を|1 ∗ · · · ∗ ∗ · · ·
···
···
··· ∗
 0 0 · · · 0 |1 ∗ · · · · · ····· · · ∗



·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
· · ∗0
A = 0······· · · 0 |1 ∗ · · · · · · ∗


···
···
···
···
· · · 0
 0
..
..
..
..
..
.. ..
.
.
.
.
.
. .

行列計算における数式処理ソフト maxima の利用について

19の形の行列に変形できる。A0 を階段行列という。A0 ですべての成分がゼロである行を除いた残りの階段状になっている部分の行数を A の階数という。これは、途中の変形の仕方によらず A によって一意に定まる数である。A の階数を rank(A) と書く。階数を用いると、連立一次方程式 (3) が解を持つための必要十分条件を簡明に rank(A) = rank(A|b)といってしまうことができる。つまり、係数行列と拡大係数行列の階数が一致するとき、そのときに限って連立一次方程式は解を持つのである。さて、このことを使って maxima
で遊ぶためのコマンドを説明しよう。コマンド echelon は、行列を階段行列に基本変形する。その構文は echelon(A); である。また、コマンド rank は行列の階数を求める。構文は rank(A); である。
まずは基本的な例から
例 18. (C1)µ A :¶ matrix([a, b], [c, d]);
a b
(D1)
c d
(C2) rank(A);(D2) 2(C3) µechelon(A);
¶
1 ab
(D3)
0 1
(C4) 
B : matrix([1,2, 3], [4, 5, 6], [7, 8, 9]);

1 2 3(D4)  4 5 6 
7 8 9(C5) rank(B);(D5) 2(C6) 
echelon(B);
1 2 3(D6)  0 1 2 
0 0 0(C7) 
C : matrix([1,2, 3], [4, 5, 6], [7, 8, 0]);
1 2 3
(D7)  4 5 6 
7 8 0(C8) rank(C);(D8) 3(C9) 
echelon(C);
1 87 0(D9)  0 1 14 
0 0 1
次に例 17 の (C9) について、可解性の判定をしてみよう。
例 19. (C1) eq1 : 3*x - 2*y = 3;
(D1) 3x − 2y = 3(C2) eq2 : x + y = 7;
(D2) y + x = 7(C3) eq3 : x - y = 9;

20

足立健朗

(D3) x − y = 9(C4) 
A : coefmatrix([eq1,eq2, eq3], [x, y]);
3 −2(D4)  1 1 
1 −1(C5) 
Ab : augcoefmatrix([eq1,eq2, eq3], [x, y]);
3 −2 −3(D5)  1 1 −7 
1 −1 −9(C6) rank(A);(D6) 2(C7) rank(Ab);(D7) 3(C8) 
echelon(A);
1 − 23(D8)  0 1 
0 0(C9) 
echelon(Ab); 
1 − 23 −1
(D9)  0 1 − 1850 01係数行列と拡大係数行列の階数が違うから当然解は存在しない。solve を使う前に調べておけば良かった！

行列計算における数式処理ソフト maxima の利用について

21

5.4. 同次連立一次方程式. さて、連立一次方程式 (4) において、b = 0 のとき、すなわち
Ax = 0 のとき、特に同次連立一次方程式と言う。このとき、rank(A) = rank(A|0) は明かであるから、同次連立一次方程式は必ず解を持つ。一般に連立一次方程式の解に含まれる任意パラメータの数は、未知数の数から係数行列の数をひいたものである。すなわち連立一次方程式 (4) の場合ならば、n − rank(A) である。したがって、n = rank(A) の場合
は任意パラメータを含まない唯一の解を持つ。このとき、特にそれが同次連立方程式ならば、唯一の解は x1 = x2 = · · · = xn = 0 であり、自明解と呼ばれる。
連立方程式 (4) において、方程式の数と未知数の数が等しい場合、つまり m = n の場合、係数行列 A は n 次正方行列になる。そして、この場合 n = rank(A) であることは det(A) 6= 0 と同値であるので、解の個数が唯一つかどうかを調べるのに A の行列式
det(A) を使うことができる。さらに、特筆すべきは、det(A) 6= 0 のとき連立一次方程式(4) の解をあからさまに
x = A−1 bと表現できることである。
例 20. (C1)2, 3], [4, 5, 6], [7, 8, 9]);
 A : matrix([1,
1 2 3(D1)  4 5 6 
7 8 9(C2) rank(A);(D2) 2(C3) solve([x+2*y+3*z=0, 4*x+5*y+6*z=0, 7*x+8*y+9*z=0], [x,y,z]);Dependent equations eliminated: (3)(D3) [[x = %R1, y = −2%R1, z = %R1]](C4) eq1 : x + 2*y + 3*z = 0;
(D4) 3z + 2y + x = 0(C5) eq2 : 4*x + 5*y + 6*z = 0;
(D5) 6z + 5y + 4x = 0(C6) eq3 : 6*x + 7*y = 0;
(D6) 7y + 6x = 0(C7) 
B : coefmatrix([eq1,eq2, eq3], [x, y, z]);
1 2 3(D7)  4 5 6 
6 7 0(C8) rank(B);(D8) 3(C9) solve([eq1, eq2, eq3], [x, y, z]);(D9) [[x = 0, y = 0, z = 0]](C10) eq4 : x + 2*y + 3*z = 10;(D10) 3z + 2y + x = 10(C11) eq5 : 4*x + 5*y + 6*z = 11;(D11) 6z + 5y + 4x = 11(C12) eq6 : 7*x + 8*y = 12;(D12) 8y + 7x = 12(C13) C : coefmatrix([eq4, eq5, eq6], [x, y, z]);

22

足立健朗




1 2 3
(D13)  4 5 6 
7 8 0(C14) Cb : augcoefmatrix([eq4,eq5, eq6], [x, y, z]);
1 2 3 −10
(D14)  4 5 6 −11 
7 8 0 −12(C15) rank(C) - rank(Cb);(D15) 0(C16) bb : 
matrix([10], [11], [12]);10(D16)  11 
12(C17) xx : C^^(-1)
. bb;

− 28
3

(D17)  29
30(C18) ££
solve([eq4, eq5, eq6],¤¤ [x, y, z]);2829x = − 3 ,y = 3 ,z = 0(D18)(C19) determinant(A);(D19) 0(C20) determinant(B);(D20) 24
(C21) determinant(C);(D21) 27
注意 4. 公式 x = A−1 b を使って解を求めるのは、solve を使うのより効率が悪い。

行列計算における数式処理ソフト maxima の利用について

23

6. 正方行列6.1. 行列式とトレース. 行列 A が n 次元正方行列であるというのは、もちろん A が (n, n)
型 行列
a11 a12 · · · a1n
 a21 a22 · · · a2n 
A=
..
.. 
...
 ...
.
. 
an1 an2 · · · ann
であることを意味する。ここで A の行列式 det(A) の定義をちゃんと書いておこう。
X
det(A) :=sgn(σ)a1σ(1) a2σ(2) · · · anσ(n)σ∈Sn

ここで Sn は n 次対称群である。maxima の determinant は正方行列の行列式を計算するものであるが、上の定義をそのまま実行するものではない。内部的には、三角行列に変形して対角成分をすべて掛け合わせたものとして計算する。これは計算量の問題による。理論的には基本変形のうち、(1),(3) のみを使って三角行列に変形すれば符合を除いては行列式は変わらないことによる。maxima においては、コマンド triangularize により、与えられた行列にたいする上三角行列を得ることができる。maxima の triangularize
は基本変形 (2) も使うので行列式の値は変形によって保存されないが、行列式がゼロの場
合は変形しても行列式ゼロである。
例 21. (C1)2, 3], [4, 5, 6], [7, 8, 9]);
 A : matrix([1,
1 2 3(D1)  4 5 6 
7 8 9(C2) 
B : triangularize(A);
1 2
3(D2)  0 −3 −6 
0 0
0(C3) determinant(A);(D3) 0(C4) determinant(B);(D4) 0(C5) 
C : matrix([1,2, 3], [4, 5, 6], [7, 8, 0]);
1 2 3
(D5)  4 5 6 
7 8 0(C6) 
D : triangularize(C);
7 8 0(D6)  0 3 42 
0 0 −27(C7) determinant(C);(D7) 27(C8) determinant(D);(D8) − 567(C9) F : matrix([a, b, c], [d, e, f], [g, h, i]);

24

足立健朗




a b c
(D9)  d e f g h i(C10) determinant(F);(D10) a (ei − f h) − b (di − f g) + c (dh − eg)(C11) G : triangularize(F);

a
b
c

af − cd(D11)  0 ae − bd00(ae − bd) i + (cd − af ) h + (bf − ce) g(C12) ¡ratsimp(
determinant(F) - col(row(G, 3), 3) );
¢
0
(D12)
ここで、col 及び row はそれぞれ行列から指定された行および列を取り出すコマンドである。col( row(M, i), j) の様に使うと、行列の M の第 (i, j) 成分を取り出すことになる。
行列式と並んで、正方行列の重要な量であるトレースを考える。定義は簡単で、すべての対角成分の総和である。maxima で行列 A のトレースを求めるには、コマンド mattrace
を用いて、mattrace(A); とタイプすればいい。使う前に load(“nchrpl”) とそれ用のパッケージを読み込む必要がある。
例 22. (C1) load("nchrpl");(D1) /usr/local/share/maxima/5.9.0/share/matrix/nchrpl.mac(C2) µA : matrix([a,b], [c, d]);
¶
a b
(D2)
c d
(C3) mattrace(A);(D3) d + a(C4) 
B : matrix([a,b, c], [d, e, f], [g, h, i]);a b c(D4)  d e f g h i(C5) mattrace(B);(D5) i + e + a(C6) ¡sum(col(row(B,i), i), i, 1, 3);
¢
i+e+a
(D6)
高校で習った Caylay-Hamilton の定理を maxima で確かめておこう。証明すべき式は、
任意の 2 次正方行列
µ
¶
a b
A=
c d
に対して、式A2 − tr(A)A + det(A)I = Oである。上の例題に引続き、
例 23. (C7)
µ I :¶ diagmatrix(2, 1);
1 0
(D7)
0 1
(C8) A^^2 - mattrace(A)*A + determinant(A)*I;

行列計算における数式処理ソフト maxima の利用についてµ−a (d + a) + ad + a2 −b (d + a) + bd + ab(D8)−c (d + a) + cd + ac d2 − d (d + a) + ad
(C9) µratsimp(D8);
¶
0 0
(D9)
0 0

25

¶

正方行列 A の行列式 det(A) とトレース tr(A) は、有用な性質をもっている。G を Aと同じ型の正則行列 (det(G) 6= 0) とするとき、det(G−1 AG) = det(A),tr(G−1 AG) = tr(A)A と B を同じ型の正方行列とするとき、det(AB) = det(A) det(B)tr(A + B) = tr(A) + tr(B)tr(AB) = tr(BA)
次数が小さいときに、maxima で確かめておこう。
例 24. (C1)µ A :¶ matrix([a, b], [c, d]);
a b
(D1)
c d
(C2) µG : matrix([s,t], [u, v]);
¶
s t
(D2)
u v
(C3) ³determinant(A) -´ ³determinant(G^^(-1)´ ³ . A . G);´³´
s(du+cs)u(bu+as)v(bv+at)t(dv+ct)(bu+as)vt(du+cs)s(dv+ct)u(bv+at)
(D3)
−
−
−
−
−
+
sv−tu
sv−tu
sv−tu
sv−tu
sv−tu
sv−tu
sv−tu
sv−tu
ad − bc(C4) ratsimp(D3);(D4) 0(C5) 
B : matrix([a,b, c], [d, e, f], [g, h, i]);a b c(D5)  d e f 
g h i(C6) 
F : matrix([o,p, q], [r, s, t], [u, v, w]);o p q(D6)  r s t 
u v w(C7) ratsimp( determinant(B) - determinant(F^^(-1) . B . F) );(D7) 0(C8) load("nchrpl");(D8) /usr/local/share/maxima/5.9.0/share/matrix/nchrpl.mac(C9) ratsimp(mattrace(A) - mattrace(G^^(-1) . A . G));
(D9) 0(C10) ratsimp(mattrace(B) - mattrace(F^^(-1) . B . F));
(D10) 0

26

足立健朗

例 25. (C1)µ A :¶ matrix([a, b], [c, d]);
a b
(D1)
c d
(C2) µB : matrix([s,t], [u, v]);
¶
s t
(D2)
u v
(C3) ratsimp(determinant(A . B) - determinant(A) * determinant(B));(D3) 0(C4) 
C : matrix([a,b,c],[d,e,f],[g,h,i]);a b c(D4)  d e f g h i(C5) 
D : matrix([o,p,q],[r,s,t],[u,v,w]);
o p q(D5)  r s t 
u v w(C6) ratsimp(determinant(C . D) - determinant(C) * determinant(D));(D6) 0例 26. (C1)µ A :¶ matrix([a, b], [c, d]);
a b
(D1)
c d
(C2) µB : matrix([s,t], [u, v]);
¶
s t
(D2)
u v
(C3) load("nchrpl");(D3) /usr/local/share/maxima/5.9.0/share/matrix/nchrpl.mac(C4) mattrace(A . B);(D4) dv + bu + ct + as(C5) mattrace(B . A);(D5) dv + bu + ct + as(C6) mattrace(A . B) - mattrace(B . A);(D6) 0(C7) 
C : matrix([a,b, c], [d, e, f], [g, h, i]);a b c(D7)  d e f 
g h i(C8) 
D : matrix([o,p, q], [r, s, t], [u, v, w]);o p q(D9)  r s t 
u v w(C10) mattrace(C . D);(D10) iw + f v + cu + ht + es + br + gq + dp + ao(C11) mattrace(D . C);(D11) iw + f v + cu + ht + es + br + gq + dp + ao(C12) mattrace(C . D) - mattrace(D . C);

行列計算における数式処理ソフト maxima の利用について

27(D12) 0
補足 1. 与えられた行列 A の行列式を求めるもう一つのコマンド newdet が用意されている。実は、newdet の方が determinant より良いアルゴリズムを使っている。determinant でうまく計算できないときには newdet を使ってみるといい。

6.2. 転置行列、対称行列、交代行列、エルミート行列、歪エルミート行列. 行列 A の行と列を入れ換えたものを A の転置行列と言って t A と書いた。maxima で、与えられた
行列 A の転置行列を求めるには、コマンド transpose を使う。
例 27. (C1)µ A :¶ matrix([a, b], [c, d]);
a b
(D1)
c d
(C2) µAt : transpose(A);
¶
a c
(D2)
b d
(C3) 
B : matrix([a,b, c], [d, e, f], [g, h, i]);a b c(D3)  d e f 
g h i(C4) 
Bt : transpose(B);a d g(D4)  b e h 
c f i(C5) µC : matrix([1,¶ 2, 3, 4], [5, 6, 7, 8]);
1 2 3 4(D5)
5 6 7 8(C6) 
Ct : transpose(C);
1 5
 2 6 (D6) 
 3 7 4 8
正方行列 A が条件 A = t A を満たすとき、A を対称行列と言う。また、A が条件
A = −t A を満たすならば A を交代行列と言う。任意の正方行列 A に対して、
A + tAA − tAB=,C=
22と置くと、B と C はそれぞれ、対称行列と交代行列になる。これを、A の対称化と交代化と言う。A = B + C のように、任意の正方行列は対称行列と交代行列の和の形に書ける。maxima には残念ながら、与えられた正方行列を対称化、交代化するコマンドは用意されていないが、maxima の持つプログラム機能を用いて、自分自身の新たなコマン
ドとして実現できる。
例 28. (C1) symm(X) := (1/2)*(X + transpose(X));(D1) symm (X) := 12 (X + TRANSPOSE (X))
(C2) µA : matrix([1,2], [3, 4]);
¶
1 2
(D2)
3 4

28

足立健朗

(C3) µsymm(A);
¶
1 52
(D3)
5
4
2
(C4) 
B : matrix([3, -1, 2], [4, -2, 1], [9, 9, 5]);3 −1 2(D4)  4 −2 1 
9 9 5(C5) 
symm(B);
3 32 112(D5)  32 −2 5 
11
5 5
2(C6) altm(X) := (1/2)*(X - transpose(X));(D6) altm (X) := 12 (X − TRANSPOSE (X))
(C7) µaltm(A);¶0 − 12(D7)
1
0
2(C8) 
altm(B);
0 − 52 − 72
(D8)  52 0 −4 
7402(C9) µA - (symm(A)+ altm(A));
¶
0 0
(D9)
0 0
(C10) B - (symm(B) + altm(B));
0 0 0(D10)  0 0 0 
0 0 0
この例では、(C1) で対称化のコマンド symm を、(C6) で交代化のコマンド altm を定義した。
与えられた行列 A に対して、各成分をその共役複素数で置き換えてできる行列を A の複素共役行列といい、A と書く。maxima で、与えられた行列の複素共役行列をもとめるには、コマンド conjugate を使う。conjugate を使うためには、その前に load(eigen);により関連パッケージをロードしておく必要がある。
例 29. (C1) load(eigen);
Warning - you are redefining the MACSYMA function EIGENVALUES
Warning - you are redefining the MACSYMA function EIGENVECTORS(D1) /usr/local/share/maxima/5.9.0/share/matrix/eigen.mac(C2) µA : matrix([1 ¶- %I, 2 + 3*%I], [%I, 4]);
1 − i 3i + 2(D2)
i4(C3) µconjugate(A); ¶i + 1 2 − 3i
(D3)
−i
4

行列計算における数式処理ソフト maxima の利用について

29

正方行列 A が条件 A = t A を満たすとき、A をエルミート行列と言う。また、A が条
件 A = −t A を満たすならば A を歪エルミート行列と言う。A を A = Re + Im と実数
部分 Re と純虚数部分に分けて考えることにしよう。A がエルミート行列ならば、Re + Im = A = t A = t Re − t Imであるから、Re は対称行列、Im は交代行列になる。また、A が歪エルミート行列ならば、Re + Im = A = −t A = −t Re + t Imであるから、Re は交代行列、Im は対称行列になる。このことから、逆に任意の正方行
列 A に対し、

Re + t Re Im − t Im+,22
Re − t Re Im + t Im
C=
+
,
2
2
と置くと、B はエルミート行列、C は歪エルミート行列になり、A = B + C が成立する。
このような、B 、C を求めるコマンドを作成してみよう。
B=

例 30. (C1) load("nchrpl");(D1) /usr/local/share/maxima/5.9.0/share/matrix/nchrpl.mac(C2) load(eigen);
Warning - you are redefining the MACSYMA function EIGENVALUES
Warning - you are redefining the MACSYMA function EIGENVECTORS(D2) /usr/local/share/maxima/5.9.0/share/matrix/eigen.mac(C3) symm(X) := (1/2)*(X + transpose(X));(D3) symm (X) := 12 (X + TRANSPOSE (X))
(C4) altm(X) := (1/2)*(X - transpose(X));(D4) altm (X) := 12 (X − TRANSPOSE (X))
(C5) ream(X) := (1/2)*(X + conjugate(X));(D5) ream (X) := 12 (X + CONJUGATE (X))
(C6) imgm(X) := (1/2)*(X - conjugate(X));(D6) imgm (X) := 12 (X − CONJUGATE (X))
(C7) herm(X) := symm(ream(X)) + altm(imgm(X));(D7) herm (X) := symm (ream (X)) + altm (imgm (X))
(C8) skew(X) := altm(ream(X)) + symm(imgm(X));(D8) skew (X) := altm (ream (X)) + symm (imgm (X))
(C9) µA : matrix([1 ¶+ %I, 2 - 3*%I], [2, 4*%I]);
i + 1 2 − 3i
(D9)
2
4i
(C10) µherm(A);¶1
2 − 3i2
(D10)
3i
+2
0
2
(C11) µskew(A); ¶i − 3i2(D11)− 3i2 4i(C12) A - (herm(A) + skew(A));

30

足立健朗

µ
(D12)

0 0
0 0

¶

注意 5. maxima で新しいコマンド (maxima、Lisp での関数のこと) を定義する最も簡単な方法はコマンドの名前 (引数１, 引数 2, . . .) := 引数１、引数 2、... の式;である。一つの式で定義できない場合は、block 文
コマンドの名前 (引数１, 引数 2, . . .) := block([局所変数の列],
引数１、引数 2、... の式１, . . . , 引数１、引数 2、... の式 N);を使う。最後の” 引数１、引数 2、... の式 N” の値がコマンドの返り値になる。

行列計算における数式処理ソフト maxima の利用について

31

7. maxima 内部における行列の扱いについて
7.1. 配列とリストと行列の関係について. リストはデータを並べてできるデータのこと、あるいはそのデータ型のことである。リストは maxima の最も基本的な構成要素であり、既存のデータ型から新しいデータ型を生成する汎用の糊の役割を果たす。maxima でのリストの表現は [a, b, ... , z] のようにデータの並びを “[“ と “]” で挟んで表す。maximaには、リストを操作するコマンドがいくつも用意されているが、詳しいことはマニュアルにゆずり、用があるときにはその都度使い方を示すことにする。
例 31. (C1) a : [1, 2, 3];(D1) [1, 2, 3](C2) a;(D2) [1, 2, 3]
配列 (array) は、いくつかの添字により指し示されたデータの集まりで、a[i, j, ..., n]の形で参照できるデータ、あるいはそのデータ型のことである。ここで、添字 i, j, ... , kの数は最大で 5 個まで可能である。添字の数が n 個の場合 n 次元配列とよばれる。配列の詳細もマニュアルにゆずることにする。
例 32. (C3) b[i] := i;(D3) bi := i(C4) b;(D4) b(C5) b[5];(D5) 5(C6) c[i,j] := i*j;(D6) ci,j := ij(C7) c;(D7) c(C8) c[3,4];(D8) 12
リストと配列は異なるデータ型であるが、maxima では、名前付きのリストの各要素にアクセスするのに、そのリストを一次元配列のように取り扱うことができる。例 31 の
a は長さ 3 のリストであるが、配列ではない。しかし、a[i] の書式で配列のようにアクセスすることができるのである。
例 33. (C9) a[1];(D9) 1(C10) a[2];(D10) 2(C11) a[3];(D11) 3
(C12) listp(a);(D13) true(C14) listp(b);(D14) false

32

足立健朗

ここで使ったコマンド listp は、その記号の値がリストなのか、そうでないのかをチェックするときに使うものである。式 listp(a); は、もし a の値がリストならば値 true を返し、もしそうでなければ値 false を返す。
リストというのは、非常に柔軟なしくみである、それ自体の入れ子を許すので、リストのリスト、リストのリストのリスト等 いくらでも考えることができる。
例 34. (C1) a : [ [1, 2, 3], [4, 5, 6], [7, 8, 9] ];(D1) [[1, 2, 3] , [4, 5, 6] , [7, 8, 9]](C2) a[1];(D2) [1, 2, 3](C3) a[2];(D3) [4, 5, 6](C4) a[3];(D4) [7, 8, 9](C5) a[1][1];(D5) 1(C6) a[1][2];(D6) 2(C7) a[1][3];(D7) 3(C8) a[2][1];(D8) 4(C9) a[2][2];(D9) 5(C10) a[2][3];(D10) 6
(C11) a[3][1];(D11) 7
(C12) a[3][2];(D12) 8(C13) a[3][3];(D13) 9(C14) a[2,3];
Wrong number of indices:[2, 3]-- an error. Quitting. To debug this try DEBUGMODE(TRUE);)
上の例で分かるように名前のついたリストのリストの要素には、それを配列の配列としてアクセスできることがわかる。上の例より a[i][j] は (a[i])[j] を意味する。しかし、2 次元配列としては扱うことができないことを示している。
行列はリストのリストでも、配列の配列でも、2 次元配列でもないが、ある行列のメン
バーにアクセスするのに、それをリストのリスト、配列の配列あるいは 2 次元配列と思ってアクセスすることが可能である。
例 35. (C1) A : matrix([a, b, c, d], [e, f, g, h], [i, j, k, l]);

行列計算における数式処理ソフト maxima の利用について
a b c d(D1)  e f g h 
i j k l(C2) A[1];(D2) [a, b, c, d](C3) A[2];(D3) [e, f, g, h](C4) A[3];(D4) [i, j, k, l](C5) A[2][3];(D5) g(C6) A[2,3];(D6) g(C7) matrixp(A);(D7) true(C8) listp(A);(D8) false(C9) matrixp(A[1]);(D9) false(C10) listp(A[1]);(D10) true(C11) B : [[a,b,c,d],[e,f,g,h],[i,j,k,l]];(D11) [[a, b, c, d] , [e, f, g, h] , [i, j, k, l]](C12) B[1];(D12) [a, b, c, d](C13) B[2];(D13) [e, f, g, h](C14) B[3];(D14) [i, j, k, l](C15) B[2][3];(D15) g(C16) B[2,3];
Wrong number of indices:[2, 3]-- an error. Quitting. To debug this try DEBUGMODE(TRUE);)(C17) matrixp(B);(D17) false(C18) listp(B);(D18) true(C19) matrixp(B[1]);(D19) false(C20) listp(B[1]);(D20) true

33

34

足立健朗

ここで使ったコマンド matrixp は、その記号の値が行列なのか、そうでないのかをチェックするときに使うものである。式 listp(a); は、もし a の値が行列ならば値 true を返し、もしそうでなければ値 false を返す。
maxima における行列型のこのような柔軟な取り扱いにより、様々な操作が可能になる。簡単な応用として、行基本変形を扱ってみよう。ここでは、まず 3 種類の n 次正則行
列 P(i,j,n), Q(a,i,n), R(a,i,j,n), (i 6= j) を定義する。(n, k) 行列 M に対して、基本変形(1) “第 i 行と第 j 行を入れ換える” ことと P(i,j,n) を M に左から掛けることは同値である。(n, k) 行列 M に対して、基本変形 (2) “第 i 行を a 倍する” ことと Q(a,i,n) を M に左から掛けることは同値である。(n, k) 行列 M に対して、基本変形 (3) “第 i 行に第 j 行の a倍を加える” ことと R(a,i,j,n) を M に左から掛けることは同値である。例で見てみよう。
まず、次のような内容のファイルをつくり、”em.mac” と名前を付けて保存する。

/* em.mac */P(i, j, n) := block([I, A],
I : ident(n),A : zeromatrix(n, n),A[i][i] : -1, A[j][j] : -1,
A[i][j] : 1, A[j][i] : 1,
I + A )$
Q(a, i, n) := block([A],A : ident(n),A[i][i] : a,A )$R(a, i, j, n) := block([I,A],
I : ident(n),A : zeromatrix(n,n),A[i][j] : a,
I + A )$
これだけ準備しておくと以下のような計算ができる。
例 36. (C1) load("./em.mac")$(C2) 
M : matrix([1,2,3],[4,5,6],[7,8,9]);

1 2 3(D2)  4 5 6 
7 8 9(C3) 
M1 : R(-4, 2, 1, 3) . M;
1 2
3

0 −3 −6 
(D3)
7 8
9
(C4) 
M2 : R(-7, 3, 1, 3) . M1;
1 2
3

0 −3 −6 (D4)
0 −6 −12行列計算における数式処理ソフト maxima の利用について

35(C5) 
M3 : R(-2, 3, 2, 3) . M2;
1 2
3(D5)  0 −3 −6 
0 0
0(C6) 
M4 : Q((-1/3),
2, 3) . M3;

1 2 3(D6)  0 1 2 
0 0 0(C7) 
echelon(M);
1 2 3(D7)  0 1 2 
0 0 0(C8) 
P(2,3,5);

1 0 0 0 0
 0 0 1 0 0 



0
1
0
0
0
(D8) 
 0 0 0 1 0 
0 0 0 0 1
(C9) P(2,4,5);

1 0 0 0 0
 0 0 0 1 0 



0
0
1
0
0
(D10) 
 0 1 0 0 0 
0 0 0 0 1
(C11) Q(3,2,5);

1 0 0 0 0
 0 3 0 0 0 



0
0
1
0
0
(D11) 
 0 0 0 1 0 
0 0 0 0 1
(C12) R(3,2,4,5);

1 0 0 0 0
 0 1 0 3 0 



0
0
1
0
0
(D12) 
 0 0 0 1 0 
0 0 0 0 1
(C13) determinant(P(2,4,5));(D13) − 1(C14) determinant(Q(3,2,5));(D14) 3
(C15) determinant(R(3,2,4,5));(D15) 17.2. 行と列の操作について. この小節では maxima で行列の行と列を扱う方法を述べよう。maxima にもともと備わっているコマンドに row 及び col がある。与えられた行列
A に対して、row(A, i); とタイプすると A の第 i 行が得られる、また、col(A, j); とすれ

36

足立健朗

ば A の第 j 列が得られる。厳密に言えば、row と col が返すものは、それぞれ (m, 1) 型
および (1, n) 型の行列である。
例 37. (C1)b, c], [d, e, f], [g, h, i]); A : matrix([a,a b c(D1)  d e f 
g h i(C2) 
col(A,
 1);a(D2)  d g(C3) matrixp(col(A, 1));
(D3) true(C4) listp(col(A, 1));
(D4) false(C5) ¡row(A, 1);
¢
a b c
(D5)
(C6) matrixp(row(A, 1));
(D6) true(C7) listp(row(A, 1));
(D7) false(C8) col(A,1)[1];(D8) [a](C9) col(A,1)[1][1];(D9) a(C10) col(A,1)[2][1];(D10) d(C11) col(A,1)[3][1];(D11) g(C12) row(A,1)[1];(D12) [a, b, c](C13) row(A,1)[1][2];(D13) b(C14) row(A,1)[1][3];(D14) c
一方、A[i]; は行列 A の i 行を返すが、これは型としては 1 次元配列であり、行列ではない。row(A, i); と A[i]; の返す結果の見ための同一性には注意が必要である。コマンド
transpose は行列の転置を返すものであるが、1 次元配列を引数として与えると (m, 1) 型の行列を返す。したがって、transpose(A[i]); は行列になる。
例 38. (C15) A[1];(D15) [a, b, c](C16) matrixp(A[1]);(D16) false(C17) listp(A[1]);(D17) true

行列計算における数式処理ソフト maxima の利用について

37(C18) transpose(A[1]);a(D18)  b c(C19) matrixp(transpose(A[1]));(D19) true
col を使わないで A の第 j 列を取り出すには、transpose(transpose(A)[j]); とタイプすればよい。
例 39. (C20) transpose(transpose(A)[1]);a(D20)  d g
以上で行、列を抽出することはできるようになった。それでは、指定の行、列に値を設定するにはどうしたらよいのか、また行基本変形についてはすでに前小節で触れたが、列基本変形についてはどうか、さらには、もっと一般的な操作、行や列の削除、追加といったことを maxima で実現するにはどうすればよいか？以下では、こういったことを考えていこう。
まず準備として、与えられた行列から指定したブロックを抽出するコマンド matblockと、指定したブロックにブロックデータを上書きするコマンド matreplace を新たに定義する。matblock(A, i1, i2, j1, j2); は、行列
a1,1 · · · a1,j1 · · · a1,j2 · · · a1,n
 ..
..
..
.. 
 .
.
.
. 


ai1,1 · · · ai1,j1 · · · ai1,j2 · · · ai1,n 
 .
..
..
.. 
A=
.
.
. 
 ..

a

·
·
·
a
·
·
·
a
·
·
·
a
i2,j1
i2,j2
i2,n 
 i2,1
 .
..
..
.. 
 ..
.
.
. 
am,1 · · · am,j1 · · · am,j2 · · · am,nから、行列ブロック

ai1,j1 · · · ...

ai1,j2
.. 
.

ai2,1

ai2,j1···を抽出するものだ。逆に、matreplace(A, B, i1, i2, j1, j2) は A の上の行列ブロック部分を、行列
b1,j2−j1+1 · · ·a1,j2−j1+1
..
..

B=
.
.

ai2−i1+1,1 · · ·ai2−i1+1,j2−j1+1

で上書きするものだ。次のような matblock と matreplace の定義を書いたファイルを用意し、“matblock.mac” と名前をつけよう。
/* matblock.mac */matblock(Mat,i1,i2,j1,j2) := block([_A, _m, _n, _i, _j],

38

足立健朗

_m : i2 - i1, _n : j2 - j1,_A : zeromatrix(_m + 1, _n + 1),
for _i : 1 thru _m + 1 do
for _j : 1 thru _n + 1 do
_A[_i][_j] : Mat[i1 + _i - 1][j1 + _j - 1],_A)$matreplace(Mat1,Mat2,i1,i2,j1,j2) := block([_A, _m, _n, _i, _j],_m : i2 - i1, _n : j2 - j1,_A : copymatrix(Mat1),
for _i : 1 thru _m + 1 do
for _j : 1 thru _n + 1 do
_A[i1 + _i - 1][j1 + _j - 1] : Mat2[_i][_j],_A)$
これらのコマンドを使うときには、事前に load(“matblock.mac”); とタイプして、定義ファイルを読み込んでおく。例をあげよう。
例 40. (C1) A : matrix([q,w,e,r,t],[y,u,i,o,p],[a,s,d,f,g],[h,i,j,k,l],[z,x,c,v,b]);
q w e r t
 y u i o p (D1) 
 a s d f g 
 h i j k l z x c v b(C2) load("./matblock.mac");(D2) ./matblock.mac(C3) 
matblock(A,2,4,2,4);u i o(D3)  s d f i j k(C4) ¡matblock(A,2,2,2,2);¢u(D4)(C5) 
matblock(A,2,5,3,4);i o
 d f 
(D5) 
 j k c v(C6) 
B : matrix([1,2,3],[4,5,6],[7,8,9]);

1 2 3(D6)  4 5 6 
7 8 9(C7) µC : matrix([1,2],[3,4]);
¶
1 2
(D7)
3 4

行列計算における数式処理ソフト maxima の利用について

39(C8) 
matreplace(A,B,2,4,2,4);
q w e r t
 y 1 2 3 p 



a
4
5
6
g
(D8) 
 h 7 8 9 l z x c v b(C9) 
matreplace(A,C,1,2,1,2);
1 2 e r t
 3 4 i o p (D9) 
 a s d f g 
 h i j k l z x c v bmatreplace を使えば、適切な大きさの行、列ベクトルを用意しておいて、それを与えられた行列に張り付けることができることがわかるであろう。注意することは、与えられた行列は matreplace で変化しないと言うことである。あくまで、そのコピーの値を変えるだけである。
列基本変形には、前小節で定義した正則行列 P(i,j,n), Q(a,i,n), R(a,i,j,n) をそのまま使うことができる。(k,n) 型行列 M に対して、(1) “第 i 列と第 1j 列を入れ換える” こととP(i,j,n) を M に右から掛けることは同値である。(2) “第 i 列を a 倍する” ことと Q(a,i,n)
を M に右から掛けることは同値である。(3) “第 i 列に第 j 列の a 倍を加える” こととR(a,j,i,n) を M に右から掛けることは同値である。例をあげよう
例 41. (C1) load("./em.mac");(D1) ./em.mac(C2) 
M : matrix([1,2,3],[4,5,6],[7,8,9]);

1 2 3(D2)  4 5 6 
7 8 9(C3) 
M1 : M . R(-2, 1, 2, 3);
1 0 3(D3)  4 −3 6 
7 −6 9(C4) 
M2 : M1 . R(-3, 1, 3, 3);
1 0
0(D4)  4 −3 −6 
7 −6 −12(C5) 
M3 : M2 . R(-2,2, 3, 3);

1 0 0
(D5)  4 −3 0 
7 −6 0(C6) 
M4 : M3 . Q((-1/3), 2, 3);
1 0 0(D6)  4 1 0 
7 2 0(C7) M . P(1, 3, 3);

40

足立健朗




3 2 1
(D7)  6 5 4 
9 8 7
行列に新たに行や列を付け加えるには、コマンド addrow と addcol を使う。addrowの書式は addrow(M, L1, L2, ...); である。ここで、M を (m, n) 型行列とするとき、L1,
L2, ... は長さ n のリスト、あるいは (l, n) 型の行列で、順に M の最下行に加えられたものが addrow の返り値となる。addcol の書式は addcol(M, L1, L2, ...); である。ここで、M を上記の通りとするとき、L1, L2 , ... は長さ m のリスト、あるいは (m, l) 型の行列
で、もしそれがリストのときは transpose(L1), transpose(L2), ... と変換されて、M の最右列に順次追加されたものが addcol の返り値となる。
例 42. (C1)
 A : matrix([a,b,c],[d,e,f],[g,h,i],[j,k,l]);a b c
 d e f (D1) 
 g h i 
j k l(C2) µB : matrix([1,2,3],[4,5,6]);
¶
1 2 3
(D2)
4 5 6
(C3) 
C : matrix([1,2],[3,4],[5,6],[7,8]);

1 2
 3 4 
(D3) 
 5 6 7 8(C4) addrow(A, C);
Incompatible structure - ADDROW//ADDCOL-- an error. Quitting. To debug this try DEBUGMODE(TRUE);)(C6) 
addrow(A, 
B);
a b c
 d e f 
 g h i (D6) 

 j k l 
 1 2 3 
4 5 6
(C7) 
addcol(A, C); a b c 1 2 d e f 3 4 (D7) 
 g h i 5 6 
j k l 7 8(C8) 
addrow(A, [1, 2, 3], [4, 5, 6]);a b c
 d e f 
 g h i (D8) 

 j k l 
 1 2 3 
4 5 6

行列計算における数式処理ソフト maxima の利用について

41(C9) 
addcol(A, [1, 3,
 5, 7], [2, 4, 6, 8]);a b c 1 2 d e f 3 4 (D9) 
 g h i 5 6 
j k l 7 8
matblock、matreplace、addrow と addcol を組み合わせて使えば、大抵の行列の切り貼りができることは理解されよう。個々の場合にどうすればいいかは読者への問題としておこう。
一つ応用を挙げておく。行列 A が
µ
¶
R S
A=
O T
の形をしているとしよう。ここで、A,R および T は正方行列とし、O の成分はすべてゼロとする。このとき、一般に

det(A) = det(R) det(S)がなり立つ。具体例で確認してみよう。
例 43. (C1)
µ R :¶ matrix([a, b], [c, d]);
a b
(D1)
c d
(C2) 
T : matrix([r,s, t], [u, v, w], [x, y, z]);r s t(D2)  u v w x y z(C3) µS : matrix([h,i, j], [k, l, m]);
¶
h i j
(D3)
k l m
(C4) 
O : zeromatrix(3,
2);

0 0
(D4)  0 0 
0 0(C5) 
A : addrow( addcol(R,S), addcol(O, T) );a b h i j
 c d k l m (D5) 
 0 0 r s t 
 0 0 u v w 
0 0 x y z(C6) determinant(R);(D6) ad − bc(C7) determinant(T);(D7) r (vz − wy) − s (uz − wx) + t (uy − vx)(C8) ratsimp(determinant(R)*determinant(T));(D8) ((ad − bc) rv + (bc − ad) su) z+((bc − ad) rw + (ad − bc) tu) y+((ad − bc) sw + (bc − ad) tv) x(C9) determinant(A);(D9) ad (r (vz − wy) − s (uz − wx) + t (uy − vx))−bc (r (vz − wy) − s (uz − wx) + t (uy − vx))

42

足立健朗

(C10) ratsimp(D8 - D9);
(D10) 0与えられた正方行列 A の (i, j) 成分 aij を含む行と列を取り除いてできる小行列 Ai,jをもとめるには、maxima に用意されたコマンド minor をつかう。書式は minor(A,i,j);である。ちなみに、行列 B = (bij ), bij := (−1)i+j det(Aji ), ∀i, j を A の余因子行列と言うのであった。A の余因子行列を求めるには、コマンド adjoint を使って、adjoint(A); とタイプすればいい。
例 44. (C1)
 A : matrix([a,b,c,d],[e,f,g,h],[i,j,k,l],[m,n,o,p]);a b c d e f g h (D1) 
 i j k l 
m n o p(C2) 
minor(A,1,1);

f g h(D2)  j k l 
n o p(C3) 
minor(A,2,3);
a b d(D3)  i j l 
m n p
例 45. (C1)µ A :¶ matrix([a, b], [c, d]);
a b
(D1)
c d
(C2) µB : zeromatrix(2,
2);
¶
0 0
(D2)
0 0
(C3) B[1][1] : (-1)^(1+1)*determinant(minor(A,1,1));(D3) d(C4) B[1][2] : (-1)^(1+2)*determinant(minor(A,2,1));(D4) − b(C5) B[2][1] : (-1)^(2+1)*determinant(minor(A,1,2));(D5) − c(C6) B[2][2] : (-1)^(2+2)*determinant(minor(A,2,2));(D6) a(C7) µ
B;
¶
d −b
(D7)
−c a
(C8) µadjoint(A);
¶
d −b
(D8)
−c a
(C9) µB - adjoint(A);
¶
0 0
(D9)
0 0

行列計算における数式処理ソフト maxima の利用について

43

さて、理論的にはA の余因子行列
A の行列式あった。maxima の adjoint と invert は、Lisp ではなく、maxima 言語を用いて書かれたユーティリティである。そのソースコードは短くて分かりやすいので、ここに挙げて
おく。
A−1 =

例 46. /* the following routines compute inverses and adjoints of matrices *//* if ratmx is false [the default] then the elements will not be converted
to cre-form */adjoint(mat):= block([adj,n], n:length(mat), adj:ident(n),if n#1 thenfor i thru n do for j thru n do
adj[i,j]:(-1)^(i+j)*determinant(minor(mat,j,i)),adj)$
invert(mat):= block([adj,ans], adj:adjoint(mat),ans:block([scalarmatrixp:true],adj/(row(mat,1).col(adj,1))),
if scalarmatrixp=true and length(mat)=1/* row(mat,1).col(adj,1) = determinant(mat) */then ans[1,1] else ans)$
複数の行と列を一挙に削除する場合には、コマンド submatrix を使う。与えられた行
列 M から 第 m1, m2, . . . 行と 第 n1, n2, . . . 列を削除した行列を得るには submatrix(m1,m2 , ..., M, n1, n2, ...); とタイプする。
例 47. (C1) A : matrix([a,b,c,d,e],[f,g,h,i,j],[k,l,m,n,p],[p,q,r,s,t],[u,v,w,x,y]);a b c d e
 f g h i j 



k
l
m
n
p
(D1) 
 p q r s t 
u v w x y(C2) 
B : submatrix(2,4,A,2,3);
a d e(D2)  k n p 
u x y

7.3. 成分の操作について. ここまでの知識から、与えられた行列 A の第 (i, j) 成分 aij を取り出すのは簡単で、A[i][j]; 、A[i,j]; を使えばいい。また、row と col を組み合わせて
使い col(row(A, i), j); とすれば、aij をただ一つの成分とする (1, 1) 型行列 (aij ) が得られる。成分だけを取り出すには col(row(A, i), j)[1][1]; あるいは col(row(A, i), j)[1,1]; とすればよい。
例 48. (C1) A : matrix([a, b, c, d], [e, f, g, h], [i, j, k, l]);

44

足立健朗




a b c d(D1)  e f g h 
i j k l(C2) A[2][3];(D2) g(C3) A[2,3];(D3) g(C4) ¡col(row(A,2), 3);¢g(D4)(C5) col(row(A, 2), 3)[1][1];(D5) g(C6) col(row(A, 2), 3)[1,1];(D6) g逆に、行列 A の第 (i, j) 成分 aij に値を設定するには、代入文 A[i][j] : a; 、A[i,j] : a;などを使えばいい。これらの方法では A そのものが変化してしまう。A そのものを変えたくないときは、A のコピーに代入文を使うか、あるいは matreplace を使って、D :matrix([a]); かつ matreplace(A, D, i, i, j, j); とすればよい。
例 49. (C1)b, c, d], [e, f, g, h], [i, j, k, l]);
 A : matrix([a,a b c d(D1)  e f g h 
i j k l(C2) A[2][3] : -3;(D2) − 3(C3) 
A;
a b c d(D3)  e f −3 h 
i j k l(C4) 
B : A;
a b c d(D4)  e f −3 h 
i j k l(C5) B[2][3] : x;(D5) x(C6) 
B;
a b c d(D6)  e f x h 
i j k l(C7) 
A;
a b c d(D7)  e f x h 
i j k l(C8) C : copymatrix(A);行列計算における数式処理ソフト maxima の利用について

45




a b c d(D8)  e f x h 
i j k l(C9) C[2][3] : g;(D9) g(C10) C;
a b c d(D10)  e f g h 
i j k l(C11) A;
a b c d(D11)  e f x h 
i j k l(C12) load("./matblock.mac");(D12) ./matblock.mac(C13) ¡D : ¢matrix([g]);g(D13)(C14) matreplace(A, D, 2, 2, 3, 3);a b c d(D14)  e f g h 
i j k l(C15) A;
a b c d(D15)  e f x h 
i j k l
注意 6. 上の例の (C4) から (D7) を見よ。手続き型言語の代入文とは違って、B : A; は B
が A と同じ値を持った A のコピーであることを意味しない。B の値が A である、あるいは A に B という別名をつけるという意味である。したがって、B の値をかえると Aの値も変化する。コピーをつくるには、コマンド copymatrix を使う必要がある。
行列の特定の成分に値を設定するもう一つの方法は、コマンド setelmx を使うことである。行列 M に対して、setelmx(a, i, j, M); とタイプすると、M の第 (i, j) 成分の値がa に設定される。このコマンドは M を変化させる。返り値は M 自身である。
例 50. (C1)
 M : matrix([a,b,c,d],[e,f,g,h],[i,j,k,l]);a b c d(D1)  e f g h 
i j k l(C2) 
setelmx(x, 2, 3, M);a b c de f x h 
(D2)
i j k l46

足立健朗

8. 線形写像と基底の変換
8.1. 線形写像、線形変換. n 個の成分がスカラー体 K に属するような列ベクトル全体の集合を Kn と書く。Kn の各ベクトルに、Kn のベクトルを 1 つづつ対応させる対応規則が与えられたとしよう。数学の言葉を使うと、集合 Kn から Km への一つの写像が与えられたとするのである。この対応規則で、x ∈ Kn が y ∈ Km に対応しているとしよう。成分
を使って、
 
 
x1
y1
.
.. 



.
x=
y=
. ,
.
xn
ym
と表すとき、y の各成分が x の成分たちの同次一次式になっているとき、すなわち対応
規則がy1 = a11 x1 + a12 x2 + · · · + a1n xny = a x + a x + · · · + a x
2
21 1
22 2
2n n

···



ym = am1 x1 + am2 x2 + · · · + amn xnの形で書けていて、しかもこの形が特定の x, y によらないとき、この対応を以前は一次写像と呼んでいた。いま、写像に f という名前をつけて、y = f (x) という記法を使うことにしよう。行列 A ∈ Mat(m, n; K)(K 上の m × n 行列全体の集合) をa11 a12 · · · a1n
 a21 a22 · · · a2n 
A=
..
.. 
...
 ...
.
. 

am1 am2 · · ·amnにより定めると、上の事は

y = f (x) = Axと表現できる。行列 A は一次写像に f によって一意的に定まり、逆に行列 A ∈ Mat(m, n; K)を任意に与えれば、上の式によって Kn から Km への一次写像 f が一つ定義されたことになる。行列 A を一次写像 f の表現行列と言う。
上の一次写像の定式化は非常に強力なものなのだが、実際の問題に適用するには不十分であることがすぐにわかる。ある問題のモデルが Kn から Km への一次写像としてとらえられるとしよう。すなわち、観測対象の集合として X と Y というものがあり、考えている問題が X から Y への写像 f によって記述されるとする。ここで、X は Kn に、
Y は Km にみたてられ、f は一次写像と考えられる。ここで問題になるのは、X あるいは Y を Kn あるいは Km に同一視する方法が一つとは限らないということである。同一視の方法がちがえば、当然 f の表現行列は違うものになるが、現象を記述している f そのものは本質的に同じものでなくてはならない。必要なことは、個々の同一視を離れた視点と、２つの特別な同一視が与えられたとき、それらの間の関係を記述する方法である。
そこで、数学の持つ抽象の力が発揮される。
まず Kn たちの持っていた性質を抽象して行くことにしよう。集合 V が K 上のベクトル空間であると言うのは、V 上に２つの演算 加法” + “ と スカラー倍 “·” が定義されていて、以下の性質が成立することを言う。
(1) v + w = w + v, ∀v, ∀w ∈ V(2) (u + v) + w = u + (v + w), ∀u, ∀v, ∀w ∈ V

行列計算における数式処理ソフト maxima の利用について

47(3) ∃1 0 : v + 0 = 0 + v = v, ∀v ∈ V(4) ∀v ∈ V, ∃1 − v ∈ V : v + (−v) = (−v) + v = 0(5) α · (β · v) = (αβ) · v, ∀α, ∀β ∈ K, ∀v ∈ V(6) (α + β) · v = α · v + β · v, ∀α, ∀β ∈ K, ∀v ∈ V(7) α · (v + w) = α · v + α · w, ∀α ∈ K, ∀v, ∀w ∈ V(8) 1 · v = v, ∀v ∈ V(9) 0 · v = 0, ∀v ∈ V
ここで、(3) の 0 を零ベクトル、(4) の −v を v の逆ベクトルと言う。(−1)v = −v であることはすぐわかる。ベクトル空間の任意の要素をベクトルと言う。以下、”·” は省略することにしよう。
抽象的なベクトル空間だけでは色々な計算ができないので、これを列ベクトル全体の集合と同一視するための言葉を与えよう。ベクトル空間 V から有限個 (n 個) の順序のついたベクトルの組 A = (a1 , a2 , . . . , an ) を選ぶことができて、これに対し任意の v ∈ V に対
して
 
v1
 v2 n
∃1 vA = 
: v = v1 a1 + v2 a2 + · · · + vn an
 ...  ∈ K
vnとなるとき、すなわち行列の積を使った記法では、
 v1v = (a1 , . . . , an ) ...  = AvA
vnと書けるとき、ベクトル空間 V は n 次元であると言い、A を V の基底と言う。また、vA
を a の A に関する成分と言う。n 次元ベクトル空間の基底の定め方は無限の可能性があるが、各基底を構成するベクトルの個数は常に n 個であることが証明できる。このあたりの詳しいことは線形代数学の教科書に譲ろう。2 通りの基底があれば、各ベクトルを成分で表す方法も 2 通りあることになるが、それらの間に成り立つ関係については後の節で解説する。
特に、 Kn 自身、n 次元ベクトル空間である。それは、標準基底 E = (e1 , e2 , . . . , en ) により、任意の列ベクトル v ∈ Kn が v1
 v2 
v=
 ...  = v1 e1 + v2 e2 + · · · + vn en = Evvnと書けることからすぐわかる。さらに、Kn においては、v = vE
が成り立っている。
すべてのベクトル空間が有限次元であるというわけではない、有限次元でないベクトル空間は無限次元であるというが、この文書では扱わない。
つぎに、一次写像を抽象化しよう。抽象ベクトル空間の間の写像を考えるので、まずは座標によらない定義が必要になる。そこで、一次写像の持っていた幾何学的、あるいは抽象代数的側面に注目しよう。幾何学的に言うと一次写像は (複素) 直線を (複素) 直線ある

48

足立健朗

いは原点に移す写像である、これを代数の言葉で表現すると次のようになる。V を K 上の n 次元ベクトル空間、 W を K 上の m 次元ベクトル空間とする。f は V から W への写像であり、性質f (αv + βw) = αf (v) + βf (w),∀α, ∀β ∈ K,∀v, ∀w ∈ V

を満たす。この性質を満たすような V から W への写像を 線形写像 と言う。線形写像という概念が一次写像を抽象化したものであることは、V と W に基底の組を入れてみれば容易に理解することができる。今、V の基底を A = (a1 , . . . , an )、W の基底をB = (b1 , . . . , bm ) としよう。線形写像 f は (f (a1 ), . . . , f (an )) ∈ W n で完全に決定される。なぜならば、任意の v ∈ V は A = (a1 , . . . , an ) たちによって、v = v1 a1 + · · · vn an = AvAと書け、f が線形写像であるから、

f (v) = v1 f (a1 ) + · · · + vn f (an )となるからである。記号の簡単のため、f (A) = (f (a1 ), . . . , f (an )) と書くことにすると、

f (v) = f (A)vAと書ける。一方、各 f (ai ) たちは W のベクトルなので、B = (b1 , . . . , bm ) によって、f (a1 ) = a11 a1 + a21 a2 · · · + am1 bmf (a ) = a a + a a · · · + a b
2
12 1
22 2
m2 m

···



f (an ) = a1n a1 + a2n a2 · · · + amn bmと書くことができる。ここで、行列 A ∈ Mat(m, n; K) をa11 a12 · · · a1n
 a21 a22 · · · a2n 
A=
..
.. 
..
 ...
.
.
. 

am1 am2 · · ·amnと定める。この行列 A は写像 f と基底の組 (A, B) によって決まる。上のことを行列の積で表現すると
f (A) = BAとなる。したがって

f (v) = BAvAを得るが、この最後の式は f (v) の B に関する成分を wB と書くことにすると、(f (v) =BwB より、)
wB = AvAを意味することがわかる。このことは、V と Kn を、 W と Km を基底対 (A, B) によって同一視したとき、線形写像 f が行列 A で定まる一次写像になっていることを意味する。
つまり、線形写像はその具体的表現では一次写像になるのである。さて、ここで出てきた
行列 A を線形写像 f の基底対 (A, B) に関する表現行列と言う。基底を取り替えたときに表現行列がどうなるかについては後の節で扱う。線形写像 f の性質は表現行列の基底の取り方に依らない性質の中に見出すことができる。

行列計算における数式処理ソフト maxima の利用について

49

8.2. 核、像、次元. ベクトル空間 V の部分集合 W がそれ自身ベクトル空間であるとき、W を V の部分ベクトル空間 (以下、単に部分空間と言うことにする) と言う。これは、
性質v, w ∈ W ⇒ αv + βw ∈ W, ∀α, ∀β ∈ Kが満たされることと言ってもよい。(注意: 一般には αv + βw ∈ V である。) 言葉で言えば、加法とスカラー倍に対して閉じた部分集合を部分空間と言うのである。V 自身と、{0} は明らかに V の部分空間であるが、これらを自明な部分空間と言う。一般に、V の任意の部分集合 S に対して、S を含むような V の部分空間のうちで、包含関係について最小のものが存在する。それを、hSi と書いて、S で生成される V の部分空間という。
集合論の記法を用いれば、
hSi = {v | v = α1 v1 + · · · + αk vk ,α1 , . . . , αk ∈ K,

v1 , . . . , vk ∈ S,

k ∈ N}

である。つまり、S の有限個の要素たちの一次結合で書ける V の要素全体の集合を hSi と書くのである。S が最初から部分空間ならば S = hSi である。ベクトル空間 V が有限次元で基底 A = (a1 , . . . , an ) を持つならば、V = h{a1 , . . . , an }i である。このときベクトル空間の次元は n であるのであるが、このことを記号で、dim V = n あるいは dimK V = nと書く。有限次元ベクトル空間 V の部分空間 W はやはり有限次元であるが、このとき
dim W ≤ dim V が成り立つ。等号成立は V = W を意味する。
この節で特に関心のあるのは、線形写像にともなう部分ベクトル空間たちである。今、V と W をともに有限次元ベクトル空間で、 dim V = n、dim W = m とする。また、
f : V → W を V から W への線形写像としよう。このとき、V の部分集合

Ker f := {v ∈ V | f (v) = 0}
を f の核 (Kernel) と言い、W の部分集合

Im f := {f (v) | v ∈ V }を f の像 (Image) と言う。核 Im f 、像 Im f は f の線形性からそれぞれ V 、 W の部分
空間になる。実際、

v, w ∈ Ker f ⇒ f (αv + βw) = αf (v) + βf (w) = 0⇒ αv + βw ∈ Ker f
及び、αf (v) + βf (w) = f (αv + βw) ∈ Im fが成り立つ。前小節で、線形写像を基底対のもとで列ベクトルの空間の間の一次写像と同一視することを説明したが、その状況の下で核と像はそれぞれ Kn と Km の部分空間と同一視される。これらを具体的に見ていこう。A = (a1 , . . . , an ) を V の基底、B = (b1 , . . . , bm )を W の基底とし、A ∈ Mat(m, n; K) を f の基底対 (A, B) に関する表現行列としよう。
v = AvA ,vA ∈ Knが v ∈ Ker f であるということは、

AvA = 0を意味する。すなわち、v の A に関する成分 vA は同次連立一次方程式 Ax = 0 の解である。この逆も真である。つまり、Ker f はその具体的表現としては、表現行列を係数行列とする同次連立一次方程式の解全体の空間 (以下、解空間という) として出現するのである。そこで、5 節の連立一次方程式の話しを全面的に適用することができる。

50

足立健朗

例 51. V = K3 , W = K2 とする。どちらも標準基底で考えよう。線形写像 f : V → W が
 
 
 
µ ¶
µ ¶
µ ¶
0
1
0
0
2
2






, f 0 =f 0 =, f 1 =
−2
−1
3
1
0
0
で与えられるものとする。このとき、f の表現行列 A は
µ
¶
2 0
2
A=
3 −1 −2である。この例では

Ker f = {v ∈ K3 |Av = 0}である。maxima を使って Ker f を調べよう。
(C1) µA : matrix([2,¶ 0, 2],[3, -1, -2]);
2 0
2
(D1)
3 −1 −2(C2) 
v : matrix([s],[t],[u]);s(D2)  t u(C3) µo :¶matrix([0],[0]);0(D3)0(C4) eq1 : (A . v)[1][1] = o[1][1];(D4) 2u + 2s = 0(C5) eq2 : (A . v)[2][1] = o[2][1];(D5) − 2u − t + 3s = 0(C6) solve([eq1, eq2],[s, t, u]);(D6) [[s = −%R1, t = −5%R1, u = %R1]]
この結果から、 * +¯
 −r ¯ −1 ¯
−5
−5r ¯ r ∈ K =Ker f =

¯


r
1 
となる。したがって、dim Ker f = 1 であることもわかる。

¤

像 Im f については、まずIm f = h{f (a1 ), . . . , f (an )}iであることに注意する。各 f (ai ) の基底 B に関する成分を wB,i と書くことにすれば、
Im f と同一視されるべき Km の部分空間は

h{wB,1 , . . . , wB,n }iで与えられることは明らかであろう。一方、

ai = AaA,i ,aA,i ∈ Kn ,i ∈ {1, . . . , n}とするとき、

aA,i = ei ,i ∈ {1, . . . , n}行列計算における数式処理ソフト maxima の利用について

51

であり、

wB,i = AaA,i = Aei ,i ∈ {1, . . . , n}となる。ここで、

fi = Aei
は行列 A の第 i 列を表すから、

A = (f1 , . . . , fn ) = (wB,1 , . . . , wB,n )が成り立つ。したがって Im f の具体的実現は

h{f1 , . . . , fn }iとなる。つまり、Im f の情報はすべて表現行列 A に含まれているわけである。線形代数学によれば、行列としての階数 rank A は A の列ベクトルたちの中で一次独立なものの最大個数に等しいことが知られているが、それは正に上の列ベクトルたちで生成される部
分空間の次元に等しい。すなわち、

dim Im f = rank Aである。この数には重要な情報があるので、線形写像 f の 階数と言って rank f で表す。

rank f := dim Im f = rank A.
例 52. 今度の例では少し抽象的にしてみよう。V は 3 次元ベクトル空間で基底 A =(a1 , a2 , a3 ) が与えられているものとする。W は 2 次元ベクトル空間で基底 B = (b1 , b2 )があてられているものとする。f は V から W への線形写像で
f (a1 ) = 3b1 + 2b2
f (a2 ) = −2b1 + b2

f (a ) = −6b − 4b312であるとしよう。このとき、f の基底対 (A, B) に関する表現行列 A はµ¶
3 −2 −6A=
2 1 −4で与えられる。f の 核、像を調べよう。
(C1) µA : matrix([3,¶ -2, 6], [2, 1, -4]);
3 −2 6(D1)
2 1 −4(C2) eq1 : 3*s - 2*t + 6*v = 0;
(D2) 6v − 2t + 3s = 0(C3) eq2 : 2*s + t - 4*v = 0;
(D3) − 4v + t + 2s = 0(C4) ££
solve([eq1, eq2], [s, t, v]);¤¤
24%R1
s = 2%R1
,
t
=
,
v
=
%R1
(D4)
7
7
(C5) rank(A);(D5) 2
これから分かるのは
¯
¾ ¿½
¾À
½
¯
2
24r
24r
2r
a1 +
a2 + ra3 ¯¯ r ∈ K =a1 +a2 + a3
Ker f =
7
7
7
7

52

足立健朗

であり、dim Ker f = 1。また、dim Im f = rank f = rank A = 2、したがって、

Im f = W¤となることである。

勝手な w ∈ W が Im f に属すかどうかというのは、次のような意味がある。w ∈ Im f
ならば、ある v ∈ V があって、w = f (v) となることである。成分を使って表すと、
wB = AvA となる。これは、非同次連立一次方程式

Ax = wBが解をもつと言うことと同値になる。一方、w ∈ Im f ということは、上述より w が
f (a1 ), . . . , f (an ) たちの一次結合で書けるということと同値である、すなわち、wB がf1 , . . . , fn の一次結合で書けることを意味する。これは、連立一次方程式 Ax = wB の係
数行列 A と 拡大係数行列 (A, wB ) の階数が同じであることと同値である。

w ∈ Im f ⇔ rank A = rank(A, wB ).
例 53. V = K3 、W = K3 としよう、簡単のためどちらも標準基底を入れて考えることにする。線形写像 f : V → W が行列


1 2 3
4 5 6 
7 8 9を使って、f (v) = Av,∀v ∈ V により定義されているとしよう。２つのベクトル
 
 
−2
−2



1 ,
1 ∈W
a=
b=
4
3

について、それぞれ Im f に属すかどうかを調べよう。
(C1) 
A : matrix([1,2, 3], [4, 5, 6], [7, 8, 9]);

1 2 3(D1)  4 5 6 
7 8 9(C2) 
a : transpose(matrix([-2,1, 4]));
−2(D2)  1 
4(C3) 
b : transpose(matrix([-2,1, 3]));
−2(D3)  1 
3(C4) rank(A);(D4) 2(C5) rank(addcol(A,a));(D5) 2(C6) rank(addcol(A,b));(D6) 3(C7) addcol(A,a);行列計算における数式処理ソフト maxima の利用について

53




1 2 3 −2
(D7)  4 5 6 1 
7 8 9 4
(C8) 
addcol(A,b); 
1 2 3 −2

4 5 6 1 (D8)
7 8 9 3
これによると、rank A = rank(A, a) より a ∈ Im f であり、rank A 6= rank(A, b) よりb∈/ Im f であることが分かる。
¤
線形写像 f : V → W の研究で最も重要な公式は次の次元公式である。

dim V = dim Ker f + dim Im f.
上述と合わせて Ker f の次元を

dim Ker f = dim V − rank f = dim V − rank Aで計算できる。Ker f は具体的表現においては 同次連立一次方程式の Ax = 0 の解空間と同一視されたから、これは解空間の次元公式でもある。
例 54. ここでも簡単のため、V = K3 , W = K3 とし、どちらにも標準基底を入れて考えよう。線形写像 f : V → W は、行列


1 2 3
A = 4 5 6 
7 8 9
により、f (v) = Av, ∀v ∈ V で定義されているものとする。
(C1) 
A : matrix([1,2, 3], [4, 5, 6], [7, 8, 9]);

1 2 3(D1)  4 5 6 
7 8 9(C2) 
v : transpose(matrix([s,t, u]));
s(D2)  t u(C3) 
AV : A . v; 
3u + 2t + s(D3)  6u + 5t + 4s 
9u + 8t + 7s
(C4) eq1 : AV[1][1];(D4) 3u + 2t + s(C5) eq2 : AV[2][1];(D5) 6u + 5t + 4s
(C6) eq3 : AV[3][1];(D6) 9u + 8t + 7s(C7) solve([eq1, eq2, eq3], [s, t, u]);Dependent equations eliminated: (3)(D7) [[s = %R1, t = −2%R1, u = %R1]]

54

足立健朗

(C8) rank(A);(D8) 2(C9) 
A;

1 2 3(D9)  4 5 6 
7 8 9(C10) load("./em.mac");(D10) ./em.mac(C11) A1 : A . R(-2, 1, 2, 3);
1 0 3(D11)  4 −3 6 
7 −6 9(C12) A2 : A1 . R(-3, 1, 3, 3);
1 0
0(D12)  4 −3 −6 
7 −6 −12(C13) A3 : A2 . R(-2,2, 3, 3);

1 0 0
(D13)  4 −3 0 
7 −6 0(C14) A4 : A3 . Q((-1/3), 2, 3);
1 0 0(D14)  4 1 0 
7 2 0
この計算から、 * +

¯
r


 1 
¯
−2
Ker f = −2r ¯¯ r ∈ R =




r
1
であり、dim Ker f = 1 であることがわかる。さらに、dim Im f = rank f = rank A = 2であり、*     + *   +
1 
2
3 
 1
 1









4 , 2
4 , 5 , 6
Im f =
=




7
0
7
8
9
であることもわかる。3 = 1 + 2 より、次元公式

dim V = dim Ker f + dim Im fの成立がたしかめられた。

¤

この節で見たように、線形写像の研究は具体的には、連立一次方程式の研究と同じである。

8.3. 基底と成分の変換. 前小節までは、基底を固定して線形写像をその表現行列から研究した。この小節では、基底を取り替えてみて、ベクトルの成分や表現行列がどのように変わるのか、そして、それらが変化しても核や像の次元が変わることはないことを見ていく。V を n 次元のベクトル空間、W を m 次元のベクトル空間とする。また f を V から
W への線形写像としよう。まず、V に 2 通りの基底 A = (a1 , . . . , an ) と A0 = (a01 , . . . , a0n )行列計算における数式処理ソフト maxima の利用について55が与えられたとしよう。基底 A0 を基準に考える。基底 A を構成する各ベクトル ai たち
も V のベクトルであるから、a01 , . . . , a0n たちの一次結合で書ける。そこで、

a1 = p11 a01 + p21 a02 + · · · + pn1 a0n a = p a0 + p a0 + · · · + p a0
2
12 1
22 2
n2 n

···



an = p1n a01 + p2n a02 + · · · + pnn a0nとしよう。ここで、行列 P ∈ Mat(n, n; K) を
p11 p12 · · · p21 p22 · · ·P := 
..
...
 ...
.
pn1 pn2 · · ·


p1n
p2n 
.. 
. 
pnn

と置くと、上のことはA = A0 Pと行列の積で表現される。 この行列 P を基底の変換 A → A0 の変換行列と言う。見方を変えると、これは恒等写像 ι : V → V の基底対 (A, A0 ) に関する表現行列と言ってもよい。rank ι = rank P = n であるから、P は正則行列である。V の任意のベクトル v に対して、その 2 通りの基底についての成分の間の関係を調べよう。

AvA = A0 vA0 (= v)より

A0 P vA = A0 vA0を得る。すなわち、

vA0 = P vA
これが、基底の変換 A → A0 にともなうベクトルの成分の変換則である。
V = Kn であるときには、基底を構成する各ベクトルも列ベクトルであるから、基底A と A0 はともに n × n 型の正則行列と見なすことができる。したがって、行列 P は、
この場合
−1
P = A0 A
により、計算することができる。
例 55. V = R3 とする。V の基底 A = (a1 , a2 , a2 ) を
 
 
2
3
a1 = −2 , a2 =  0  ,1−1とする。また、別の V の基底 A0 = (a01 , a02 , a02 ) を
 
 
1
3
0
0



a1 = 1 , a2 = 1 ,
0
3

 
1
a1 = 1
2
1
a01 = −2
4

56

足立健朗

としよう。このとき、行列として


2
3 1
A = −2 0 1 ,1 −1 2
1 3 1
A0 = 1 1 −2
0 3 4

であり、基底の変換 A → A0 の変換行列 P は P = (A0 )−1 A で与えられる。
(C1) 
A1 : matrix([2, 3, 1], [-2, 0, 1], [1, -1, 2]);2
3 1
(D1)  −2 0 1 
1 −1 2(C2) 
A2 : matrix([1,3, 1], [1, 1, -2], [0, 3, 4]);

1 3 1
(D2)  1 1 −2 
0 3 4
(C3) 
P : A2^^(-1) . A1;
31
37 −13(D3)  −13 −15 6 
10
11 −4(C4) 
A1 - (A2 . P);
0 0 0

0 0 0 (D4)
0 0 0
上の計算から31
37 −13P = −13 −15 6 
10
11 −4であることがわかる。ここで、

−1
v= 3 6の基底 A 及び A0 それぞれについての成分を調べて、基底の変換にともなう成分の変換則の成立をチェックしておこう。v の成分 vA と vA0 はそれぞれ、

vA = A−1 v,で計算できる。
(C5) 
v : transpose(matrix([-1,3, 6]));

−1

3 
(D5)
6
(C6) 
v1 : A1^^(-1)
. v;

4
− 19
(D6)  − 20 
19
49
19

(C7) v2 : A2^^(-1) . v;

−1vA0 = A0 v

行列計算における数式処理ソフト maxima の利用について

57




−79
(D7)  34 
−24(C8) 
P . v1;
−79

34 
(D8)
−24
(C9) 
v2 -P . v1;
0(D9)  0 0
上の計算より、

 4
− 19

,
vA = − 20
19
49
19




−79
vA0 =  34 
−24

であり、

vA0 = P vA¤が成立していることが確認された。

W にも 2 通りの基底 B と B0 をいれてみよう。基底の変換 B → B0 の変換行列を
Q ∈ Mat(n, n; K) とする。もちろん、B = B0 Q
であり、Q は正則行列である。線形写像 f : V → W の 2 通りの基底対 (A, B)、(A0 , B0 )についての表現行列を、それぞれ、A、A0 としよう。このとき、

f (A) = BA,

f (A0 ) = B0 A0である。f の線形性を使うと、

B0 QA = BA
= f (A)= f (A0 P )= f (A0 )P= B0 A0 P
すなわち、

QA = A0 P
あるいは、

A = Q−1 A0 P,
A0 = QAP −1を得る。これが、表現行列の変換則である。前に正方行列のところで言及したように、行列の階数は正則行列を掛けても変わらないことから、rank A = rank A0が結論される。これにより、核や像の次元は、どんな基底対で計算しても同じになることが保証される。

58

足立健朗

例 56. 例 55 の状況に、さらに、W = K2 を考える。W にも２つの基底 B = (b1 , b2 ) とB0 = (b0 1 , b0 2 ) を与えよう。ただし、
µ ¶
µ ¶
µ ¶
µ ¶
3
1
2
−1
0
0
b1 =
, b2 =, b1 =, b2 =21
1−2とする。このとき、B と B0 は行列として、
µ
¶
µ
¶
−1 3
1 2
0
B=
, B =
2 1
1 −2
と書け、基底の変換 B → B0 の変換行列 Q は Q = B0 −1 B で与えられる。今、線形写像
f : V → W の基底対 (A, B) に関する表現行列がµ¶
2 0 −1F =4 1 −5で与えられているとして、f の基底対 (A0 , B0 ) に関する表現行列 F 0 を求めてみよう。
(C1) 
A1 : matrix([2, 3, 1], [-2, 0, 1], [1, -1, 2]);2
3 1
(D1)  −2 0 1 
1 −1 2(C2) 
A2 :matrix([1, 3, 1], [1, 1, -2], [0, 3, 4]);
1 3 1
(D2)  1 1 −2 
0 3 4
(C3) 
P : A2^^(-1) . A1;
31
37 −13(D3)  −13 −15 6 
10
11 −4(C4) µB1 : matrix([-1,3], [2, 1]);
¶
−1 3
(D4)
2 1
(C5) µB2 : matrix([1,2], [1, -2]);
¶
1 2
(D5)
1 −2
(C6) µQ : B2^^(-1)
. B1;
¶
1
2
2
(D6)
− 34 12(C7) µ
F1 : matrix([2,0, -1], [4, 1, -5]);¶
2 0 −1(D7)
4 1 −5(C8) µ
F2 : Q . F1 . (P^^(-1));
¶
41
− 223
− 495
38
38
19
(D8)
− 45
− 181
− 23
76
76
19
(C9) µQ^^(-1) . F2
¶ . P;
2 0 −1(D9)
4 1 −5行列計算における数式処理ソフト maxima の利用について

上の計算より、

µ 1
Q=

及び、
0

F = QF P
を得る。

−1

2

− 34

2

59

¶

1
2

µ 223
¶
41
− 38 − 495
38
19
=
45
− 181
− 23
− 76
76
19
¤

60

足立健朗

9. ベクトルの内積
9.1. ベクトルの内積について. ベクトルの内積とは何だろうか？ 高校で習った知識によると、それは 2 つの平面あるいは空間ベクトル a と b に対して、一つの実数 a • b が決
まり、これを a と b の内積と言ったのであった。そして、その定義は θ を a と b の間の角とするとき、a • b := |a||b| cos θ
で与えられるものだった。あるいは、例えば、成分で a = (a1 , a2 , a3 )、b = (b1 , b2 , b3 ) と与えられているならば、a • b = a1 b1 + a2 b2 + a3 c3としてよいことも習ったであろう。さて、内積を一般のベクトル空間でも定義するためには、すこし反省してみることが必要である。勝手なベクトル空間におけるベクトルの長さ
とか、2 つのベクトルの間の角とか、そんなものがあるのかどうかわからないし、成分にしたって何を基準に決めてよいかわからないからである。
そこで、平面あるいは空間ベクトルの内積の持つ形式的な性質を抽象して、一般のベクトル空間における内積を定義することにしよう。与えられたベクトル空間 V 上の内積とは、直積空間 V × V から、スカラー体 K への写像 g であって、任意の a1 , a2 , b1 , b2 ∈ Vと、任意のスカラー α1 , α2 , β1 , β2 ∈ K に対して、

g(α1 a1 +α2 a2 , β1 b1 +β2 b2 ) = α1 β1 g(a1 , b1 )+α1 β2 g(a1 , b2 )+α2 β1 g(a2 , b1 )+α2 β2 g(a2 , b2 )及び、任意の a, b ∈ V に対して

g(a, b) = g(b, a)を満たすものである。一般には、さらに次の性質を附化して考えることが普通である。

g(a, a) ≥ 0
であり、等号成立は a = 0 のときに限る。この最後の性質を満たす内積を特に正定値内積というが、単に内積と言ってしまうことが多い。それで、普通困らないのであるが、例
えば相対性理論などでは、正定値でない内積を扱うのでそういう場合は意識的にこの性質を除いて考える必要がある。
さて、この文書では maxima で実際に計算できる話しに注意を向けているので、考えるベクトル空間はすべて有限次元とし、スカラー体 K は複素数体 C または 実数体 R とする。ベクトル空間 V が n 次元であるとすると、V の基底 A = {a1 , a2 , · · · , an } を一つ決めることができる。前節のように、任意のベクトル a ∈ V は、この基底を用いて

a = α1 a1 + α2 a2 + · · · + αn anと書け、このとき


α1
 α2 

aA = 
 ...  ∈ K
αn行列計算における数式処理ソフト maxima の利用について

61

を a の A に関する成分という。V の任意のベクトル a、b ∈ V を考え、それらの A に関する成分がそれぞれ
 
 
α1
β1
 α2 
 β2 

aA = 
bA = 
 ...  ,
 ... αnβnであるとする。このとき、g(a, b) は内積の性質により、
g(a, b) = α1 β1 g(a1 , a1 ) + α1 β2 g(a1 , a2 ) + · · · + αn βn g(an , an )
n X
n
X
=
αi βj g(ai , aj )i=1 j=1となる。ここで、簡単のため

gij := g(ai , aj ),∀i, ∀j ∈ {1, 2, . . . , n}と書くことにする。この gi,j たちを行列のかたちに並べて
g11 g12 · · · g1n
 g21 g22 · · · g2n 
G := 
.. . .
.. 
 ...
.
.
. 

gn1 gn2 · · ·gnnと書くことにする。内積 g の性質は、基底 A を定めれば、この行列 G によって完全に決まってしまう。G を g の A に関する計量行列という。g を G を使って表現すると次のようになる。
 
g11 g12 · · · g1n
β1 g21 g22 · · · g2n   β2 .g(a, b) = (α1 , α2 , . . . , αn ) .. . .
.. 
 ...
.
.
.   .. 

gn1 gn2 · · ·

gnn

βn

= t aA GbA .
ところで、g(a, b) = g(b, a), ∀a, ∀b ∈ V であるから、

gij = gji ,∀i, ∀j ∈ {1, 2, . . . , n}が成り立つ。すなわち、

G = tG
であり、K = C ならば G はエルミート行列、K = R ならば G は対称行列である。
話しを少し簡単にして、V = Kn の場合を考えよう。V の基底として、標準基底 E =(e1 , e2 , . . . , en ) を使うことにすれば、任意の a ∈ V に対して、a = aE となるから、ベク
トルとその E に関する成分を区別しないで用いることにする。任意に n 次エルミート行
列 G を与えるとg(a, b) := t aGb, ∀a, ∀b ∈ V
により、(正定値とは限らない) 内積 g : V × V −→ K が一つ定まる。このことを使って、与えられた計量行列 G から内積 g を作る maxima のコマンドを作ってみよう。以下の内容を “metrics.mac” と名付けてファイルに保存する。

62

足立健朗

/* metrics.mac */load(eigen)$geom(a, b, G) := block([ta],ta : conjugate(transpose(a)),
ratsimp(ta . G . b))$
コマンド geom に引数として与える行列 G はエルミート行列あるいは対称行列でないと内積にはならない。早速テストしてみよう。
例 57. (C1) load("metrics.mac");
Warning - you are redefining the MACSYMA function EIGENVALUES
Warning - you are redefining the MACSYMA function EIGENVECTORS(D1) metrics.mac(C2) µG : matrix([1,¶ 2 - %i], [2 + %i, 3]);
1
2−i
(D2)
i+2
3
(C3) µa :¶matrix([s], [t]);s(D4)t(C5) µb : ¶matrix([x], [y]);x(D5)y(C6) geom(a, b, G);(D6) (3t + (2 − i) s) y + ((i + 2) t + s) x
(C7) geom(b, a, G);(D7) (3t + (i + 2) s) y + ((2 − i) t + s) x
(C8) geom(a, a, G);(D8) 3t2 + 4st + s2
(C9) geom(b, b, G);(D9) 3y 2 + 4xy + x2
(C10) kill(a,b);(D10) DONE
(C11) µa :¶matrix([1], [2]);
1
(D11)
2
(C12) µb : matrix([2- 3 * %i], [3 + %i]);¶2 − 3i(D12)
i+3(C13) geom(a, b, G);(D13) 41 − 6i(C14) geom(b, a, G);(D14) 6i + 41
(C15) geom(a, a, G);(D15) 21

行列計算における数式処理ソフト maxima の利用について

63(C16) geom(b, b, G);(D16) 77
さて、計量行列 G が単位行列 I である場合:

1 0 0 ···
0 1 0 · · ·0 0 1
G=I=
..
 ... ...
.

0 0 00 0 0 ···
すなわち、(
gij = δij =


0 0
0 0

0 0
,
.. 
.

1 0
0 1

1 (i = j)
0 (i =
6 j)

であるとき、これに対応する内積 g を Kn の標準内積という。特に K = C のとき Cn の標準内積をエルミート内積、K = R のとき Rn の標準内積をユークリッド内積ともいう。
maxima には、標準内積を計算するコマンド innerproduct が用意されている。また省略形 inprod も使える。
例 58. (C1) load("metrics.mac")
Warning - you are redefining the MACSYMA function EIGENVALUES
Warning - you are redefining the MACSYMA function EIGENVECTORS(D1) metrics.mac(C2) 
a : transpose(matrix([s,t, u]));
s(D2)  t u(C3) 
b : transpose(matrix([x,y, z]));
x(D5)  y z(C6) inprod(a, b);(D6) uz + ty + sx(C7) 
I : ident(3);
1 0 0(D7)  0 1 0 
0 0 1(C8) geom(a, b, I);(D8) uz + ty + sx(C9) kill(a, b)(D9) DONE(C10) ¡a : matrix([s,t, u]);
¢
s t u
(D10)
(C11) ¡b : matrix([x,y, z]);
¢
x y z
(D11)

64

足立健朗

(C12) inprod(a, b);(D12) uz + ty + sx(C13) geom(a, b, I);(D13) uz + ty + sx(C14) kill(a, b);(D14) DONE
(C15) a : [s, t, u];(D18) [s, t, u](C19) b : [x, y, z];(D16) [x, y, z](C17) inprod(a, b);(D19) uz + ty + sx(C20) geom(a, b, I);(D20) uz + ty + sx
inprod を使う前には “load(eigen)” でライブラリを読み込んでおく必要がある。標準内
積は正定値内積である。標準内積 g に対し、値 g(a, b) を
ha, bi := g(a, b)と表すことにする。Kn の任意のベクトル a に対して、a の絶対値 (大きさ、長さ) ||a|| を||a||2 := ha, aiにより定義する。また、K = R のとき、Rn の任意のゼロベクトルでない２つのベクトルa, b に対して、a と b の間の角 θ (0 ≤ θ ≤ π) を

cos θ :=

ha, bi||a|| ||b||で定義する。
以上、標準計量をもとにして、ベクトルの絶対値、間の角を定義した。話しの順序は逆であるが、これで標準内積が高校で学んだ平面と空間ベクトルの内積の一般化になってい
ることがわかった。
一般のベクトル空間 V の話しに戻る。V と V 上の内積 g が与えられれば、標準内積の場合のアナロジーとして、||a||2 := g(a, b),cos θ :=g(a, b)||a|| ||b||により、g に関するベクトルの絶対値、2 つのベクトルの間の角を定義できることに気付くだろう。V = Kn の場合に、与えられた計量行列 G から決まる内積 g に関するベクトルの絶対値を求めるコマンド absolute と、2 つのベクトルの間の角の cosine を求めるコ
マンド anglcos を定義しよう。以下のプログラムを “metrics.mac” に追加しよう。

absolute(a, G) := ratsimp(sqrt(geom(a, a, G)))$
anglcos(a, b, G) := block([denom],
denom : absolute(a, G) * absolute(a, G),

行列計算における数式処理ソフト maxima の利用について

ratsimp(geom(a, b, G) / denom))$
これを、使った計算例を見てみよう。
例 59. (C1) load("metrics.mac");
Warning - you are redefining the MACSYMA function EIGENVALUES
Warning - you are redefining the MACSYMA function EIGENVECTORS(D1) metrics.mac(C2) 
I : ident(4);
1 0 0 0
 0 1 0 0 (D2) 
 0 0 1 0 
0 0 0 1
(C3) len(a) := absolute(a, I);(D3) len (a) := absolute (a, I)(C4) ancos(a, b) := anglcos(a, b, I);(D4) ancos (a, b) := anglcos (a, b, I)(C5) 
a : transpose(matrix([1,-2, 0, -1]));1
 −2 

(D5)  0 −1(C6) 
b : transpose(matrix([2,1, 3, -1]));

2
 1 

(D6)  3 −1(C7) √
len(a);(D7) 6(C8) √
len(b);(D8) 15(C9) ancos(a, b);1(D9) √6√
15(C10) inprod(a, b)/ (len(a) * len(b));1(D10) √6√
15
(C11) kill(a,b);(D11) DONE
(C12) a : transpose(matrix([s,t, u, v]));

s
 t 

(D12)  u v(C13) b : transpose(matrix([w, x, y, z]));

65

66

足立健朗




w
 x 

(D13)  y z(C14) √
len(a);(D14) v 2 + u2 + t2 + s2
(C15) plen(b);(D15)
z 2 + y 2 + x2 + w 2(C16) inprod(a, b);(D16) vz + uy + tx + sw
(C17) ancos(a, b);
vz+uy+tx+sw√(D17) √
v 2 +u2 +t2 +s2

z 2 +y 2 +x2 +w2

注意しておきたいことは、内積も絶対値も角も、ベクトル空間 V の特別な基底の取り方に依存しない概念であるということである。(Kn の標準計量ですら、標準基底によらずに定まるものである。) 基底の取り方を変えたときに変わるのは計量行列の成分だけである。V に 2 通りの基底 A と B を与える。2 つのベクトル a と b の A に関する成分をそ
れぞれ aA 、bA とし、 a と b の B に関する成分をそれぞれ aB 、bB とする。また、内積 g の A に関する計量行列を GA 、B に関する計量行列を GB とする。内積 g(a, b) そのものは基底の取り方に依らないのであるから、ta

A GA bA =ta G b
B B B(= g(a, b))が成立しているはずである。今、基底の変換 A → B の変換行列を P とすると、

A = BPであり、

aB = P aA ,bB = P bAである。それゆえ、ta

B GB bB =ta tP G P b = ta G b
A
B
A
A A A

を得る。a および b は任意であるから、
tP G

B P = GAあるいは、
t P −1 G P −1 = G
ABを得る。これが、基底の変換にともなう計量行列の成分の変換則である。
ベクトル空間 V に内積 g が与えられているとき、V と g をセットにして考え、(V, g) を内積空間あるいは、計量ベクトル空間という。特に、(Cn , h, i) をエルミート空間、(Rn , h, i)
をユークリッド空間という。
正値でない内積の例も扱っておこう。Rn の標準基底 E に関して、計量行列


1
O


−1

M =
..


.

O

−1

行列計算における数式処理ソフト maxima の利用について

67

により定まる Rn の内積 η を (1, n − 1) 型のローレンツ内積という。特に、R4 のとき (1, 3)型のローレンツ内積をミンコフスキー内積という。そして、内積空間 (R4 , η) を R1,3 と書いてミンコフスキー空間という。これは、アインシュタインの特殊相対性理論を数学的に定式化するときに出てくるものである。任意のベクトル a ∈ R1,3 に対して、η に関する自己内積 η(a, a) はもはや非負とは限らない。η(a, a) < 0 のとき、a を空間的ベクトル、
η(a, a) = 0 のとき、a を光的ベクトル、そして、η(a, a) > 0 のとき、a を時間的ベクトルという。普通の意味では絶対値を定義することはできないので、かわりにミンコフスキー
ノルムというものを考える (記号は絶対値と同じ記号を使う)。 その定義は、
p −η(a, a) 空間的な場合||a|| := 0
光的な場合pη(a, a)
時間的な場合である。ベクトル a の第 1 成分が正のとき、ベクトル a は未来的であると言うことにしよう。2 つの時間的ベクトル a, b ∈ R1,3 に対し、b − a が未来的時間的ベクトルであるとき、a から b をみた双曲角 µ を

cosh µ :=η(a, b)||a|| ||b||で定義する。
maxima でこれらを計算させるコマンドを作ってみよう。MinkMatrix はミンコフ
スキー内積の計量行列をあたえる。mink はミンコフスキー内積、tminknorm、sminknorm はそれぞれ時間的ベクトル、空間的ベクトルのミンコフスキーノルムを計算する。
minknorm は引数のベクトルが数値として具体的に与えられている場合に、そのミンコ
フスキーノルムを計算する。ただし、文字式をいれるとエラーになる。２つの時間的ベク
トルの間の双曲角の hyperbolic cosine を計算するには hypangcos を用いる。これらを、以下のようにコーディングしよう。

MinkMatrix : ematrix(4, 4, 2, 1, 1) - ident(4)$mink(a, b) := geom(a, b, MinkMatrix)$tminknorm(a) := sqrt(mink(a,a))$sminknorm(a) := sqrt(-mink(a,a))$minknorm(a) := block(if mink(a, a) < 0 thensqrt(-mink(a, a))elsesqrt(mink(a, a)))$hypangcos(a, b) := ratsimp(mink(a, b) / (tminknorm(a) * tminknorm(b)))$
再び、”metrics.mac” に追加しておこう。計算例は以下の通りである。

68

足立健朗

例 60. (C1) load("metrics.mac");
Warning - you are redefining the MACSYMA function EIGENVALUES
Warning - you are redefining the MACSYMA function EIGENVECTORS(D1) metrics.mac(C2) 
MinkMatrix;

1 0
0
0
 0 −1 0
0 (D2) 
 0 0 −1 0 
0 0
0 −1(C3) 
a : transpose(matrix([s,t, u, v]));

s
 t 

(D3)  u v(C4) 
b : transpose(matrix([w,x, y, z]));

w
 x 

(D4)  y z(C5) mink(a, b);(D5) − vz − uy − tx + sw
(C6) √
tminknorm(a);(D6) −v 2 − u2 − t2 + s2
(C7) √
sminknorm(a);(D7) v 2 + u2 + t2 − s2(C8) hypangcos(a, b);
vz+uy+tx−sw√(D8) − √−v 2 −u2 −t2 +s2−z 2 −y 2 −x2 +w2(C9) minknorm(a);MACSYMA was unable to evaluate the predicate:
ERREXP 1# 0: minknorm(a=MATRIX([s],[t],[u],[v]))(metrics.mac line 10)-- an error. Quitting. To debug this try DEBUGMODE(TRUE);)(C10) kill(a, b);(D10) DONE
(C11) a : transpose(matrix([5,1, 0, -2]));

5
 1 

(D11)  0 
−2(C12) b : transpose(matrix([6,2, -1, 2]));

6
 2 

(D12)  −1 
2

行列計算における数式処理ソフト maxima の利用について

69(C13) mink(a, b);(D13) 32
(C14) mink(a, a);(D14) 20
(C15) mink(b, b);(D15) 27
(C16) minknorm(a);√(D16) 2 5(C17) minknorm(b);√(D17) 3 3(C18) hypangcos(a, b);√(D18) 3√16
3 5(C19) c : transpose(matrix([1,2, -1, 3]));

1
 2 

(D19)  −1 
3(C20) mink(c, c);(D20) − 13(C21) √
sminknorm(c);(D21) 13
(C22) √
minknorm(c);(D22) 139.2. Gram-Schmidt の直交化. 内積空間 (V, g) において、与えられた有限個のベクトルの順序列 (a1 , . . . , ak ) に適当な操作を施して、直交系 (a01 , . . . , a0k ) を作る方法が、GramSchmidt の直交化法と言う名で知られている。ここで、 (a01 , . . . , a0k ) が直交系であると言うのは、(
g(a0i , a0i ) 6= 0 ∀i ∈ {1, . . . , k}g(a0i , a0j ) = 0 ∀i 6= ∀j ∈ {1, . . . , k}が成り立つことを意味する。条件 g(v, w) = 0 が成り立つとき、v と w は (g に関して) 直交すると言い v ⊥ w と表すこともある。(a01 , . . . , a0k ) が直交系となるためには、(a1 , . . . , ak ) も一定の条件を満たす必要がある。その条件は

h{a1 , . . . , ak }i = h{a01 , . . . , a0k }iである。これは、適当な操作が構成ベクトルの一つに他の構成ベクトルのスカラー倍を加えるという、可逆な操作を使うことによる要請である。ここで注意したいのは、{a01 , . . . , a0k }は線形独立になることで、これは次のことが成り立つことを意味する。

α1 a01 + · · · + αk a0k = 0
であるとする。∀i ∈ {1, . . . , k} について、上の式の両辺と a0i との内積を取ってみると、
内積 g の性質より、
αi g(a0i , a0i ) = 0

70

足立健朗

を得るが、g(a0i , a0i ) 6= 0 であったから、αi = 0, ∀i でなければならない。このことは、さらに (a01 , . . . , a0k ) が h{a01 , . . . , a0k }i の基底であることも示している。個数が等しいので、{a1 , . . . , ak } は線形独立でなければならない。以上から、議論の出発点を正確に述べることができるようになった。
改めて、(V, g) を n 次元内積空間とし、(a1 , . . . , ak ) (k ≤ n) を V の k 個の線形独立なベク
トルよりなる順序列とする。ただし、各ベクトルはヌルベクトルでは無い (g(ai , ai ) 6= 0, ∀i)とする。これに対し、順序列 (a01 , . . . , a0k ) を以下のように帰納的に構成する。まず、

a01 := a1とする。(a01 , . . . , a0l ) (l ≤ k − 1) が直交系であるように構成できたとして、al+1 0 を

al+1 0 := al+1 −g(al+1 , a01 ) 0g(al+1 , a0l ) 0
a
−
·
·
·
−
a
g(a01 , a01 ) 1g(a0l , a0l ) lとする。このとき、∀i ∈ {1, . . . , l} に対して、

g(al+1 0 , a0i ) = g(al+1 , a0i ) −g(al+1 , a0i )g(a0i , a0i )00g(ai , ai )=0を得る。また、{a01 , . . . , a0l , al+1 } の一次独立性より、

g(al+1 0 , al+1 0 ) 6= 0を得る。したがって、(a01 , . . . , al+1 0 ) はやはり直交系になる。この操作を繰り返し、最終的に、直交系 (a01 , . . . , a0k ) を得る。以上のようにして、(a1 , . . . , ak ) から (a01 , . . . , a0k ) を作る方法を Gram-Schmidt の直交化法と言う。
さらに、内積 g が正定値であるときには、

a00i := p1a0i ,g(a0i , a0i )∀i ∈ {1, . . . , k}とすれば、(a001 , . . . , a00k ) は、いわゆる正規直交系と呼ばれるものになる。すなわち、(
g(a00i , a00i ) = 1 ∀i ∈ {1, . . . , k}g(a00i , a00j ) = 0 ∀i 6= ∀j ∈ {1, . . . , k}を満たす。(a1 , . . . , ak ) から (a001 , . . . , a00k ) を求めるこの方法のことを Gram-Schmidt の正規直交化法と言う。
さて、大学初年級においては実用上は Rn あるいは Cn に標準内積を入れたもの、すな
わちユークリッド空間あるいはエルミート空間の場合に限定して考察すれば十分である。
以下では、標準内積の場合に限定して、Gram-Schmidt 法を maxima で実行する手順について説明しよう。
maxima のライブラリコマンド gramschmidt は、引数として、リストのリストをとり、値としてやはりリストのリストを返す。すなわち、各ベクトルをリストとして表現
し、それらの並びでベクトルの順序列を表すのである。maxima においては、行列はリストのリストとしてもあつかえるので、引数に行列を与えることも可能である。この場
合、行列の各行ベクトルが並んだ順序列と解釈される。gramschmidt を使うには、事前に “load(eigeni);” とタイプして、ライブラリを読み込んでおかねばならない。

行列計算における数式処理ソフト maxima の利用について

71

例 61. (C1) load(eigeni);
Warning - you are redefining the MACSYMA function EIGENVALUES
Warning - you are redefining the MACSYMA function EIGENVECTORS(D1) /usr/local/share/maxima/5.9.0/share/matrix/eigeni.mac(C2) £gramschmidt([[2,3],[1,-2]]);£¤¤27(D2) [2, 3] , 37
,
−
13
13
(C3) a1 : [1, 2, 3];(D3) [1, 2, 3](C4) a2 : [4, 5, 6];(D4) [4, 5, 6](C5) a3 : [7, 8, -1];
(D5) [7, 8, −1](C6) hgramschmidt([a1,
a3]);
h 2
ia2,
£ 5 25 5 ¤i
2 3 3
23(D6) [1, 2, 3] , 7 , 7 , − 7 , − 3 , 3 , − 3(C7) 
A : matrix([1, 2, 3], [4, 5, 6], [7, 8, -1]);
1 2 3(D7)  4 5 6 
7 8 −1(C8) hgramschmidt(A);
h 2
i £
¤i
5 25
5
(D8) [1, 2, 3] , 273 , 37 , − 23
,
−
,
,
−
7
3 3
3
ベクトルの順序列を、列ベクトルの並びとして、あるいはそれを行列として表現したもので与えたい場合がよく起こる。あるいは、正規直交化した結果が欲しい場合もあるだろう。maxima の標準ライブラリコマンドには直接そのような取り扱いのできるコマンドがないので、新たなコマンドとして自作することにしよう。以下のプログラムを “gson.mac”と言うファイルとして保存しよう。

/* gson.mac */load(eigeni)$gso(A) := block([m,n,i,B,C,D],m : length(A),B : transpose(A),n : length(B),C : gramschmidt(B),D : zeromatrix(n, m),
for i:1 thru n do
D[i] : C[i],transpose(D))$
gson(A) := block([m,n,i,B,C,D],m : length(A),

72

足立健朗

B : transpose(A),n : length(B),C : gramschmidt(B),D : zeromatrix(n, m),
for i:1 thru n do
D[i] : unitvector(C[i]),transpose(D))$
ここで定義したコマンド gso と gson はともに行列を引数として、行列を値として返す。gso は引数の行列を列ベクトルの順序列と解釈して Gram-Schmidt 法を適用したものを返す。gson は Gram-Schmidt の正規直交化法を適用したものを返す。コマンド定義に使われている maxima のコマンド unitvector は引数にリストを取り、それをベクトルと思って正規化 (単位ベクトル化) したものを値として返すものである。
例 62. (C1) load("gson.mac");
Warning - you are redefining the MACSYMA function EIGENVALUES
Warning - you are redefining the MACSYMA function EIGENVECTORS(D1) gson.mac(C2) µA : transpose(matrix([2,-1], [-1, 2]));
¶
2 −1
(D2)
−1 2
(C3) µratsimp(gso(A));
¶
2 35
(D3)
−1 65
(C4) Ã
ratsimp(gson(A));
!
√

(D4)

√2
5
− √15

5
5
√
2 5
5

(C5) ha[1] : transpose(D4)[1];i(D5) √25 , − √15(C6) ha[2] : itranspose(D4)[2];
√
√
5 2 5
(D6)
,
5
5
(C7) inprod(a[1], a[2]);(D7) 0(C8) 
A : transpose(matrix([1,2, 3], [4, 5, 6], [7, 8, -1]));
1 4 7(D8)  2 5 8 
3 6 −1(C9) 
ratsimp(gso(A));

5
−
1 12
7
3
10 
(D9)  2 37
3
3 − 67 − 53(C10) ratsimp(gson(A));

行列計算における数式処理ソフト maxima の利用について(D10) 

√1
14
√2
14
√3
14

√
4√7
√7 √3
3 7
21 √
√
2 3 7
− 21

√



√ √
2 3
3√
− 3√32− 3√32

73(C11) B : transpose(matrix([1,0, -2, 1], [2, 3, 3, -1], [0, 1, -1, -1]));

1
2
0
 0
3
1 

(D11)  −2 3 −1 
1 −1 −1(C12) rank(B);(D12) 3(C13) ratsimp(gso(B));
17
50
1
− 113
6
80
 0

3
113 
(D13)  −2 4 − 90 
3113
1 − 16 − 130
113(C14) ratsimp(gson(B));
√
√

17
5 113
√1
√ 6
√
−
6 √
113
113 3
√
 6

3 6
8 113
√
√
 0

113
113

√ √3
√
(D14)  − √23 3 113 
4 6
√
−


113
6
3 √
113
√
13 √
113
√1
√6
−
−
6
6 113
113 3
引数の行列が正則行列である場合、コマンド gson の返り値は、次の小節で議論する直交行列、あるいはユニタリ行列になる。A = (a1 , . . . , an ) が正則行列であるとする。このとき、A00 = (a001 , . . . , a00n ) = gson(A) を正規直交化された行列とすると、
 
t a00
1
t A0 =  .. 
 . 
t a00
n

であり、


t A0 A0 =t a00 a00
1 2
t a00 a00
2 2

..
.

···
···
...

t a00 a00
1 n
t a00 a00 
2 n

t a00 a00
n 1

t a00 a00
n 2···

t a00 a00
n n




..
.

となる。ここで、∀i, ∀j ∈ {1, . . . , n} に対して、(
t a00 a00 = ha00 , a00 i =jii jであることに注意すれば、を得る。したがって、

t a00 a00
1 1
t a00 a00
 2 1

..
.




1 (i = j)
0 (i =
6 j)

t A00 A00 = I
nt A00 = (A00 )−1であり、K = R ならば A00 は直交行列、K = C ならば A00 はユニタリ行列である。

74

足立健朗

例 63. (C1) load("gson.mac");
Warning - you are redefining the MACSYMA function EIGENVALUES
Warning - you are redefining the MACSYMA function EIGENVECTORS(D1) gson.mac(C2) µA : matrix([1,2], [3, 4]);
¶
1 2
(D2)
3 4
(C3) Ã
AA : ratsimp(gson(A));
!
√

(D3)

√1
10
√3
10

3√5
5 √2
− 5√52

√1
10
√
3√5
5 2

√3
10
√
− 5√52

10
10
√
3√2
2 5

3 10
10
√
− 2√25(C4) Ã
transpose(AA);!
(D4)(C5) Ã
AA^^(-1);
√
√
(D5)

!

(C6) µratsimp(AA^^(-1)- transpose(AA));
¶
0 0
(D6)
0 0
(C7) µratsimp(transpose(AA). AA);
¶
1 0
(D7)
0 1
(C8) 
B : matrix([1+%i, 2, -%i], [0, 2-%i, 1-2*%i], [0, 0, 2+3*%i]);i+1
2−i
2 − i 1 − 2i (D8)  000
3i + 2(C9) rank(B);(D9) 3(C10) BB : ratsimp(gson(B));

i+1
√
0
0
 2
√
0 
(D10)  0 − i−2

5
3i+2
√
0
0
13
(C11) conjugate(transpose(BB));

1−i
√
0
0
 2
√
0 
(D11)  0 − −i−2

5
2−3i
√
0
0
13
(C12) BB^^(-1);
√

2
0
0
√
 i+1

(D12)  0 − i−25
0 
√
13
0
0
3i+2
(C13) ratsimp(conjugate(transpose(BB)) - BB^^(-1));

行列計算における数式処理ソフト maxima の利用について

75




0 0 0
(D13)  0 0 0 
0 0 0(C14) ratsimp(conjugate(transpose(BB)). BB);

1 0 0(D14)  0 1 0 
0 0 1(C15) µC : matrix([1+%i,
2-%i], [3+2*%i, 4*%i]);¶
i+1 2−i(D15)
2i + 3 4i(C16) rank(C);(D16) 2(C17) ÃCC : ratsimp(gson(C));√√ !
11 5i−10
5
i+1
√
√
−
15
√ 5 51√(D17)
5 5i−3 5
2i+3
√
15

√
5 51

(C18) µratsimp(conjugate(transpose(CC))- CC^^(-1));
¶
0 0
(D18)
0 0
(C19) µratsimp(conjugate(transpose(CC)). CC);
¶
1 0
(D19)
0 1

9.3. 直交行列とユニタリ行列. ベクトル空間 V から V 自身への線形写像のことを特にV 上の線形変換と言う。V に内積 g が与えられているとしよう。すなわち、(V, g) は内積空間とする。V 上の線形変換の中でも、内積を保つものは興味深い研究対象になる。
そこで、内積を保つ線形変換という概念を正確に記述しておこう。改めて (V, g) を内積空間とし、f : V → V を V 上の線形変換とする。f が内積 g を保つとは、条件
g(f (a), f (b)) = g(a, b),∀a, ∀b ∈ Vが成り立つことを言う。
f が内積 g を保つ線形変換であるとき、任意の V の基底 A に関して、f の基底対
(A, A) に関する表現行列 F と g の A に関する計量行列 G の間には何らかの関係があるのだろうか。以下では、それを考えてみる。V の任意のベクトル a と b の基底 A に関する成分をそれぞれ aA , bA ∈ Kn としよう。ただし、n = dim V である。このとき、f がg を保つという条件を A に関する成分表示で表すと、
t (F aA )G(AbA ) =t a (t F GA)b = t a GbAAAAとなる。aA および bA は任意であるから、欲しかった条件は簡潔に
t F GF = Gと表現できる。この関係式を満たす行列 F を計量行列 G に関する直交行列と言う。
ここで、特別な場合について考えておこう。V = Kn とし、内積として標準内積 h•, •iを考える。基底として標準基底 E を考えよう。標準内積の標準基底に関する計量行列は

76

足立健朗

単位行列 In であった。そこで、標準内積を保つ任意の線形変換 f の基底対 (E, E) に関する表現行列 F は In に関する直交行列となることが分かる。すなわち、
tF I

nF =

tF F = Inが成り立つ。この条件を満たす行列は特に重要なので特別に名前が与えられている。K = Rのとき、t
F F = Inを満たす行列 F を単に直交行列と言う。K = C のとき、tF F = Inを満たす行列 F をユニタリ行列と言う。また、これに対応して Kn 上の標準内積を保つ
線形変換 f を K = R のとき直交変換、K = C のときユニタリ変換と言う。
再び、一般の正定値内積空間 (V, g) で話しをする。V の基底 A = (a1 , . . . , an ) が正規直交基底であるということを、条件(
1 (i = j)g(ai , aj ) =0 (i 6= j)が満たされることと定義する。V の別の基底 B = (b1 , . . . bn ) もまた正規直交基底であるとしよう。このとき、任意のベクトル v ∈ V の A、B それぞれに関する成分 vA 、vBを考える。前節で述べたように、基底の変換 A → B の変換行列を P とすれば、

vB = P vAがなりたつ。このとき、行列 P は係数体 K が R ならば直交行列であり、C ならばユニタリ行列である。なぜならば、

hei , ej i = g(ai , aj ) = g(bi , bj ),∀i, ∀j

であり、それゆえ

hP vA , P wA i = hvB , wB i = g(v, w) = hvA , wA iが任意の v, w ∈ V に対して成立するから、上述のように
tP P = Inを得るからである。
例 64. (C1) load("gson.mac");
Warning - you are redefining the MACSYMA function EIGENVALUES
Warning - you are redefining the MACSYMA function EIGENVECTORS(D1) gson.mac(C2) A : ratsimp(gson(matrix([1,2,3],[4,5,6],[7,8,0])));
√
√
 1

3 11
√
√2
−
11
66
2 3
√
√


11
√2
(D2)  √466

11
3
√
√
11
2
7
√
− 11 − 2√3
66(C3) B : ratsimp(gson(matrix([1,0,1],[0,1,2],[1,3,-1])));

行列計算における数式処理ソフト maxima の利用について



√1
2


(D3)  0√2− 23√11

√1
2

√
√2
11
√
3√ 2
2 11

√

11
11
√
3 11
11
√
− 1111(C4) 
P : √A^^(-1)√ √. B; √

(D4) 




66
4 √
33 2
2
√ √
2 √11
− 33

2 66
13 √
66 √
11
− 5112
√
2√ 3
3 11

√1
3
√
− i−1
3

3i−3
√ 3
3 5 √
√
√ 3
− 2 33i−
5

77

√
11 66
121
7
√
√11
2 3√ 11
11 2(C5) 
ratsimp(transpose(P)
. P);

1 0 0
(D5)  0 1 0 
0 0 1(C6) 
ratsimp(transpose(A)
. A);

1 0 0
(D6)  0 1 0 
0 0 1(C7) 
ratsimp(transpose(B)
. B);

1 0 0
(D7)  0 1 0 
0 0 1(C8) Ã
C : ratsimp(gson(matrix([1+%i,2], [2-%i, 3+2*%i])));
√
√ !
11
7i−13
7
i+1
√
− 7√587√√(D8)
10 7i+4
7
√
√
− i−2
7
7 58
(C9) Ã
D : ratsimp(gson(matrix([1,%i], [1-%i, 2-%i])));
!
√
√
(D9)

(C10) Q : ratsimp(C. √D);
√ √
√
√
√
√ √ √
√ √
√ √ √
√ √ 
3 5 58+9 3 5( 3 58−24 √3)√i+ 3 58+2 3(2 3 5 58+37 3 √5)i+4√−
3 7 58√ √58 √ √

√ √ √
√ √
√
√
√ √15 7 √√ √(D10) 3 58−14 35 3 5 58+2 3 5)i−5 3 5 58+24 3 5( 3 58−6 3√)i−2
(
√
√ √
−
3 7 5815 7 58(C11) µratsimp(conjugate(transpose(Q))
. Q);
¶
1 0
(D11)
0 1
(C12) µratsimp(conjugate(transpose(C))
. C);
¶
1 0
(D12)
0 1
(C13) µratsimp(conjugate(transpose(D))
. D);
¶
1 0
(D13)
0 1
この例では A, B, P が直交行列、C, D, Q がユニタリ行列である。

¤

ユニタリ行列あるいは直交行列の行列式については次が成り立つ。A をユニタリ行列とすると、|t A||A| = |t AA| = |In | = 1.

78

足立健朗

したがって、|A| は絶対値が 1 の複素数であることがわかる。また、A が直交行列ならば、|A|2 = |t A||A| = |t AA| = |In | = 1
より、|A| = ±1 であることがわかる。行列式が 1 であるユニタリ行列、直交行列を、それぞれ特殊ユニタリ行列、特殊直交行列と言う。
例 65. (C1) load("gson.mac");
Warning - you are redefining the MACSYMA function EIGENVALUES
Warning - you are redefining the MACSYMA function EIGENVECTORS(D1) gson.mac(C2) A : ratsimp(gson(matrix([1,2,3],[4,5,6],[7,8,0])));
√
√
 1

3 11
√
√2
−
11
66
2
3
√
√


11
√2
(D2)  √466

11
3
√
√
11
2
7
√
− 11 − 2√3
66(C3) ratsimp(determinant(A)^2);(D3) 1(C4) B : ratsimp(gson(matrix([1,%i,2-%i],[3+%i,
4*%i, 2],[-1+%i, 2-3*%i, %i])));
√
√
√
√
 1

67 193i−125
13i−
13
193
√
√
√
− 13 193−
13√√√ 193 113√ √
11 13i+1713
5 193i−38
193
i+3
√
√
(D4)  13

13 193 √
193
113√
√
√
24 13i−39
15 193i−
13
193
i−1
√
√
√
− 13 193− 193 113
13
(C5) detb : ratsimp(determinant(B));7i−8(D5) − √113(C6) ratsimp(conjugate(detb)*detb);(D6) 1

2 次の正方行列が直交行列になる条件を考えよう。
µ
¶
a b
A=
c d
が直交行列であるとすると、t AA = I2 より、
µ
¶µ
¶ µ 2
¶ µ
¶
a c
a b
a + c2 ab + cd
1 0
=
=
b d
c d
ba + dc b2 + d2
0 1を得る。この条件を読みとくと、µ¶
cos θ − sin θA=
sin θ cos θµまたは

cos η sin η
sin η − cos η¶のいずれかの形に書けることがわかる。それぞれが、行列式 1 および −1 に対応している。幾何学的には前者は平面の角 θ の回転を表す。後者は 平面の反転をあらわす行列
µ
¶
0 1
1 0
と前者のタイプの行列の積として表すことができる。

行列計算における数式処理ソフト maxima の利用について

79

例 66. (C1)-sin(th)], [sin(th), cos(th)]);¶µ A : matrix([cos(th),cos th − sin th(D1)
sin th cos th(C2) µtranspose(A) . A;¶
sin2 th + cos2 th0(D2)0
sin2 th + cos2 th(C3) µtrigsimp(D2);
¶
1 0
(D3)
0 1
(C4) trigsimp(determinant(A));(D4) 1(C5) µB : matrix([cos(phi),sin(phi)], [sin(phi), -cos(phi)]);¶cos ϕ sin ϕ(D5)
sin ϕ − cos ϕ(C6) µtranspose(B) . B;¶
sin2 ϕ + cos2 ϕ0(D6)0
sin2 ϕ + cos2 ϕ(C7) µtrigsimp(D6);
¶
1 0
(D7)
0 1
(C8) trigsimp(determinant(B));(D8) − 1(C9) µT : matrix([0,1], [1, 0]);
¶
0 1
(D9)
1 0
(C10) µC : B . T;¶
sin ϕ cos ϕ(D10)− cos ϕ sin ϕ(C11) µtrigsimp(transpose(C)
. C);
¶
1 0
(D11)
0 1
(C12) trigsimp(determinant(C));(D12) 1
上の計算例で使った maxima のコマンド trigsimp は、三角関数を含む式の表現を簡約化するコマンドである。
¤
正定値でない内積を持つ内積空間の例として、Rn にローレンツ内積 η を入れたもの、
すなわちローレンツ空間 R1,n−1 = (Rn , η) を考えよう。(n = 4 のときはミンコフスキー空間。) ローレンツ内積 η の標準基底 E に関する計量行列 J は、



1


J =


0
−1

0

..

.






−1

80

足立健朗

であった。R1,n−1 上の線形変換 f でローレンツ内積 η を保つものを考える。そのような
f の基底対 (E, E) に関する表現行列を F とすると、先に述べた議論から F は条件t

F JF = Jを満足する。一般にこの条件を満たす行列をローレンツ行列と言う。基底対 (E, E) に関する表現行列がローレンツ行列であるような R1,n−1 上の線形変換をローレンツ変換と言う。
R1,n−1 の基底 A = (a1 , . . . , an ) が条件(i = j = 1)
1
η(ai , aj ) = −1 (i = j 6= 1)0(i 6= j)を満たすとき、A を R1,n−1 のローレンツ基底と言う。これは、正定値内積空間の正規直交基底に対応する概念である。いま R1,n−1 の 2 通りの基底 A = (a1 , . . . , an )、B = (b1 , . . . , bn )がともにローレンツ基底であるとしよう。このとき、基底の変換 A → B にともなう変換
行列を P とすれば、直交行列のときと同様の議論で、条件t

P JP = Jが成り立つことがわかる。すなわち P はローレンツ行列になる。この基底の変換にともなう座標の変換公式
vB = P vA
をローレンツ座標変換と言う。F を任意のローレンツ行列としよう。定義より、|F |2 |J| = |t F ||J||F | = |t F JF | = |J|であり、|J| 6= 0 であるから、|F |2 = 1
すなわち、|F | = ±1 を得る。
2 次のローレンツ行列を考えよう。

µ
¶
a b
A=
c d

がローレンツ行列であるとすると、定義から
µ
¶µ
¶µ
¶ µ
¶
a c
1 0
a b
1 0
=
b d
0 −1
c d
0 −1
が成り立つ。すなわち、

µ 2
¶ µ
¶
a − c2 ab − cd
1 0=
ab − cd b2 − d20 −1が成り立つ。これより、A はµ¶µ¶cosh µ sinh µcosh µ −sinh µ,,sinh µ cosh µsinh µ − cosh µµ¶− cosh µ −sinh µ,sinh µcosh µのいずれかの形に書けることがわかる。v = tanh µ とすると、v
1, sinh µ = √cosh µ = √1 − v2
1 − v2µ− cosh µ sinh µsinh µ − cosh µ¶行列計算における数式処理ソフト maxima の利用についてであるから、v

µ

√ 1
1−v 2
√ v
1−v 2

81

(|v| < 1) を任意パラメータとして、
√ v
1−v 2
√ 1
1−v 2

¶

µ
,

√ 1
1−v 2
√ v
1−v 2

√ −v
1−v 2
√ −1
1−v 2

¶

µ
,

√ −1
1−v 2
√ v
1−v 2

√ −v
1−v 2
√ 1
1−v 2

¶

µ
,

√ −1
1−v 2
√ v
1−v 2

√ v
1−v 2
√ −1
1−v 2

¶

のいずれかの形と言ってもいい。
例 67. (C1) load("metrics.mac");
Warning - you are redefining the MACSYMA function EIGENVALUES
Warning - you are redefining the MACSYMA function EIGENVECTORS(D1) metrics.mac(C2) µA : matrix([cosh(mu),sinh(mu)], [sinh(mu), cosh(mu)]);¶cosh µ sinhµ(D2)
sinhµ cosh µ(C3) µJ : matrix([1,0], [0, -1]);
¶
1 0
(D3)
0 −1
(C4) µtranspose(A) . J . A;¶
cosh2 µ − sinh2 µ0(D4)0sinh2 µ − cosh2 µ(C5) µtrigsimp(D4);
¶
1 0
(D5)
0 −1
(C6) µv0 :¶transpose(matrix([t, x]));t(D6)x(C7) µv1 : A . v0;¶
sinhµx + cosh µt(D7)cosh µx + sinhµt(C8) geom(v0, v0, J);(D8) t2 − x2(C9) ¡geom(v1, v1, J);¢¡¢(D9) sinh2 µ − cosh2 µ x2 + cosh2 µ − sinh2 µ t2
(C10) trigsimp(D9);(D10) t2 − x2(C11) µB : matrix([1/sqrt(1-v^2),
v/sqrt(1-v^2)],[v/sqrt(1-v^2), 1/sqrt(1-v^2)]);
¶

(D11)

√ 1
1−v 2
√ v
1−v 2

√ v
1−v 2
√ 1
1−v 2

(C12) Ãtranspose(B) . J . B; !
v2
1
0
2 − 1−v 2
1−v
(D12)
v2
1
0
− 1−v
2
1−v 2
(C13) µratsimp(D12);
¶
1 0
(D13)
0 −1
(C14) v2 : B . v0;

82

足立健朗

µ
(D14)

t
√ vx
+ √1−v
2
1−v 2
√ x
√ tv
+
1−v 2
1−v 2

¶

(C15) geom(v2 , v2, J);(D15) t2 − x2
行列 F がユニタリ行列、(直交行列、ローレンツ行列) ならば、その逆行列 F −1 = t F もユニタリ行列 (直交行列、ローレンツ行列) である。また、F, G がともに同じ型のユニタリ
行列 (直交行列、ローレンツ行列) ならば、積 F G もまたユニタリ行列 (直交行列、ローレンツ行列) である。このことを、n 次ユニタリ行列 (直交行列、ローレンツ行列) 全体の集合は群をなすと言う。n 次ユニタリ行列全体のなす群を U (n) と書いて、n 次ユニタリ群と言う。同様に、n 次の直交行列、ローレンツ行列全体のなす群をそれぞれ O(n)、O(1, n − 1)と書いて n 次直交群、ローレンツ群と言う。特に、特殊ユニタリ行列全体、および特殊直交行列全体の集合も群をなすが、それらをそれぞれ SU (n)、SO(n) と書いて、n 次特殊ユニタリ群、特殊直交群と言う。

行列計算における数式処理ソフト maxima の利用について

83

10. 固有値と固有ベクトル10.1. 固有値と固有方程式. V を K 上のベクトル空間、f を V 上の線形変換とする。Vのあるベクトル v(6= 0) に対して、あるスカラー λ ∈ K が存在して、
f (v) = λvとなるとき、λ を f の固有値と言い、また v を λ に属する固有ベクトルと言う。
V の部分集合 Vλ を
Vλ := {v ∈ V | f (v) = λv}により定義する。このとき、

v, w ∈ Vλ , α, β ∈ K⇒ f (αv + βw) = αf (v) + βf (w) = α(λv) + β(λw) = λ(αv + βw)⇒ αv + βw ∈ Vλとなるから、Vλ は V の部分ベクトル空間である。V を λ に属する固有空間と言う。
V は有限次元 (n = dim V ) であるとし、基底 A = (a1 , . . . , an ) が与えられているとする。このとき、f (v) = λv という条件式を、基底 A に関する成分を使って表現しよう。
それは、f の基底対 (A, A) に関する表現行列を A とすると、

AvA = λvAと書ける。ここで、vA ∈ Kn は v の A に関する成分である。

hλi := {a ∈ Kn | Aa = λa}とすると、これは Vλ に 1 対 1 に対応する Kn の部分ベクトル空間であり、Vλ を基底 Aに関する成分で表現したものになっている。hλi を λ に属する Kn の固有空間と言う。
hλi 6= {0} (すなわち、Vλ 6= {0}) という条件は、λ に属する固有ベクトルが少なくとも一つ存在するということを意味し、 が固有値であるということを意味する。これは、a ∈ Kn についての同次連立一次方程式(A − λIn )a = 0が非自明解を持つという条件に他ならない。同次連立一次方程式の理論から、この方程式が非自明解をもつための必要十分条件は rank(A − λIn ) < n、あるいは、行列式の言葉で、(5)|A − λIn | = 0と表現される。逆に、これを λ に関する条件式と見ると、λ がこの方程式の解ならば、λは f の固有値であり、λ に属する固有ベクトルが存在することになる。方程式 (5) を行列
A の固有方程式、あるいは、線形変換 f の (基底 A に関する) 固有方程式と言う。ここ
で、次のことを注意しておく。f の固有値は f のみによって決まり、基底 A、したがって表現行列 A によらない。
与えられた正方行列 A に対して、その固有方程式を求める maxima のコマンドは
charpoly である。行列 A の変数 x に関する固有方程式を “charpoly(A, x);” とタイプして求めることができる。
例 68. (C1)
-1, 0], [3, 1, -2], [0, 1, 1]);
 A : matrix([2,
2 −1 0(D1)  3 1 −2 
0 1
1

84

足立健朗

(C2) ¡charpoly(A ¢, x);(D2) (1 − x)2 + 2 (2 − x) + 3 (1 − x)(C3) expand(D2);(D3) − x3 + 4x2 − 10x + 9
固有値を求めるには、もちろん、charpoly で求めた固有方程式を solve で解いてもよいのであるが、maxima にはそれ用の専用コマンド eigenvalues が用意されている。
eigenvalues は、正方行列を引数として取り、リストのリストを値として返す。返り値の最初の部分リストは固有値のリストで、次の部分リストは各固有値に対応する固有方程式の解としての重複度のリストである。
例 69. (C1)
1, -1], [1, 1, 0], [-1, 0, 1]);
 A : matrix([2,2 1 −1(D1)  1 1 0 −1 0 1(C2) eq1 : charpoly(A, x);(D2) 2x + (1 − x)2 (2 − x) − 2(C3) solve([eq1], [x]);(D3) [x = 3, x = 1, x = 0]
(C4) eigenvalues(A);Warning - you are redefining the MACSYMA function EIGENVALUES
Warning - you are redefining the MACSYMA function EIGENVECTORS(D4) [[3, 1, 0] , [1, 1, 1]]
この例では、3, 1, 0 が行列 A の固有値であり、それぞれの重複度が 1, 1, 1 であることを示している。
¤
固有値が求まったら、次は各固有値に属する固有ベクトルを求めることになるが、それはもちろん、同次連立一次方程式(A − λIn )a = 0を解くことにより求まる。その解は固有空間 hλi の次元に相当する個数の任意パラメー
タを含む形でもとまる。
例 70. 例 69 の行列 A の固有値 3, 1, 0 に属する固有ベクトルをそれぞれ求めてみよう。
(C1) A : matrix([2, 1, -1], [1, 1, 0], [-1, 0, 1]);2 1 −1
1 1 0 (D1)
−1 0 1(C2) 
xx : transpose(matrix([x, y, z]));
x

y 
(D2)
z
(C3) 
eq3 : (A - 3*ident(3)). xx;−z + y − x(D3)  x − 2y 
−2z − x(C4) solve([eq3[1][1], eq3[2][1], eq3[3][1]], [x, y, z]);Dependent equations eliminated: (3)行列計算における数式処理ソフト maxima の利用について

85(D4) [[x = −2%R1, y = −%R1, z = %R1]](C5) 
eq1 : (A - 1*ident(3)). xx;−z + y + x
x(D5) 
−x(C6) solve([eq1[1][1], eq1[2][1], eq1[3][1]], [x, y, z]);Dependent equations eliminated: (3)(D6) [[x = 0, y = %R2, z = %R2]](C7) 
eq0 : A . xx; −z + y + 2xy+x(D7) 
z−x(C8) solve([eq0[1][1], eq0[2][1], eq0[3][1]], [x, y, z]);Dependent equations eliminated: (3)(D8) [[x = %R3, y = −%R3, z = %R3]](C9) 
xx3 : transpose(matrix([-2*r,-r, r]));
−2r(D9)  −r r(C10) A . 
xx3 - 3 * xx3;0(D10)  0 
0(C11) xx1 
: transpose(matrix([0, s, s]));0(D11)  s s(C12) A . 
xx1 - 1 * xx1;
0(D12)  0 
0(C13) xx0 : transpose(matrix([t, -t, t]));t(D13)  −t t(C14) A . 
xx0 - 0 * xx0;0(D14)  0 
0
実は maxima には、直接固有ベクトルを求めてくれるコマンドが用意されている。コ
マンド eigenvactors は正方行列を引数に取り、リストのリストを値として返す。返り値の最初の部分リストは eigenvalues の返り値と同じもので、固有値のリストとその重複度のリストからなるリストのりすとである。第 2 以降の部分リストが各固有値に対応する固有ベクトルを表す。固有ベクトルは任意パラメータに適当な値を代入して得られる数ベク

86

足立健朗

トルとして返される。eigenvectors を使用するには事前に “load(eigen)” によりライブラ
リを読み込んでおく必要がある。
例 71. (C1) load(eigen)
Warning - you are redefining the MACSYMA function EIGENVALUES
Warning - you are redefining the MACSYMA function EIGENVECTORS(D1) /usr/local/share/maxima/5.9.0/share/matrix/eigen.mac(C2) 
A : matrix([2, 1, -1], [1, 1, 0], [-1, 0, 1]);2 1 −1(D2)  1 1 0 −1 0 1(C3) £eigenvectors(A);£¤¤(D3) [[3, 1, 0] , [1, 1, 1]] , 1, 12 , − 21 , [0, 1, 1] , [1, −1, 1](C4) eigenvalues(A);(D4) [[3, 1, 0] , [1, 1, 1]]
ここで、最後の式は
 
 
 
1
0
1
 1  ,  1  , −1 ,
2
−1
1
− 21
が、それぞれ、固有値 3, 1, 0 に対応する固有ベクトルであることを示している。

¤

目的によっては、eigenvector の返り値であるリスト形式ではなくて、別の形式で答が欲しい場合があるだろう。特に、次小節で使うために、固有ベクトルを列ベクトルで表し、それを行列の形に並べたものを返すようなコマンドが欲しい。そのために、新しいコ
マンドを自作することにしよう。次のプログラムを “gson.mac” に追加しよう。

eigmat(A) := block([i,n,B,C],n : length(A),B : eigenvectors(A),C : zeromatrix(n, n),for i:1 thru n do
C[i] : B[i+1],transpose(C))$
コマンド eigmat は引数の正方行列の固有ベクトルを eigenvectors で求めたものを列ベクトルにして、行列に並べたものを返す。
例 72. (C1) load("gson.mac");
Warning - you are redefining the MACSYMA function EIGENVALUES
Warning - you are redefining the MACSYMA function EIGENVECTORS(D1) gson.mac(C2) 
A : matrix([2, 1, -1], [1, 1, 0], [-1, 0, 1]);2 1 −1(D2)  1 1 0 −1 0 1

行列計算における数式処理ソフト maxima の利用について

87(C3) 
eigmat(A); 
1 0 1(D3)  12 1 −1 
− 12 1 1
一般に、n 次正方行列の固有方程式は n 次の多項式であるから、n ≥ 5 ならば当然解けない可能性が大きい。

10.2. 対称行列の対角化. 与えられた正方行列 A に対して、適当な正則行列 P を選んで、


0
λ1






λ1
−1

P AP = 
.


..






0
λn
の形にできるとき、行列 A は P により対角化可能であると言う。
上の状況が成り立っているとき、P = (p1 , p2 , . . . , pn ) を P の列ベクトルへの分解とす
れば、上の式の両辺に P を左から掛けて、



0
λ1






λ1


AP = P 
...







0
λn
を得る。すなわち、(Ap1 , Ap2 , . . . , Apn ) = (λ1 p1 , λ2 p2 , . . . , λn pn )であるが、これは、各 i ∈ {1, . . . , n} に対して、

Api = λi piが成立することを意味する。したがって、λ1 , . . . , λn はすべて A の固有値であり、p1 , . . . , pn
はそれぞれの固有値に属する A の固有ベクトルとなる。P が正則であると言う条件から、(p1 , . . . , pn ) は Kn の基底である。
V を K 上の n 次元ベクトル空間とし、f : V → V を V 上の線形変換とする。V の一つの基底 A = (a1 , . . . , an ) が与えられており、f の基底対 (A, A) に対する表現行列がA であるとしよう。いま、上のように A はある正則行列 P により対角化可能であるとする。基底 A による V と Kn の同一視は、各 ai ∈ V に ei ∈ Kn を対応させることにより得られている。この対応の逆対応により、各ベクトル pi ∈ Kn に対応するベクトル vi ∈ V を考えよう。(p1 , . . . , pn ) が Kn の基底であったから、V = (v1 , . . . , vn ) は、それに対応する V の基底となる。V における基底の変換 A → V は、Kn における基底

88

足立健朗

の変換 E = (e1 , . . . , en ) → P = (p1 , . . . , pn ) に対応しており、その変換行列はともに、
P −1 E = P −1 In = P −1 で与えられることがわかう。このことから、f の基底対 (V, V) に関する表現行列は、変換公式にしたがい、



0
λ1






λ1
−1
−1 −1
−1

P A(P ) = P AP = 
.


..






0
λn
で与えられることがわかる。つまり、線形変換 f のある基底対 (A, A) に関する表現行列
A が対角化可能であるということは、f の固有ベクトル v1 , . . . , vn から構成される基底V が存在するということを意味する。逆に、そのような基底 V = (v1 , . . . , vn ) が存在す
れば、
f (vi ) = λi vi , ∀i ∈ {1, . . . , n}より、 f の 基底対 (V, V) に関する表現行列は




0
λ1






λ1


.


.
.






0
λn
で与えられることがわかる。
そこで、f が対角化可能であるということを、固有ベクトルからなる基底が存在することと定義する。これは、特別な基底の取り方にはよらない、f のみに依存する性質である。
さて、(V, g) を R 上の正定値内積空間としよう。V 上の線形変換 f : V → V が対称変換であるとは、条件
g(f (v), w) = g(v, f (w)), ∀v, ∀w ∈ V
成り立つことである。今、A = (a1 , . . . , an ) を V の正規直交基底とし、 A を f の 基底
対 (A, A) に関する表現行列としよう。すると上式より、

g(f (ai ), aj ) = g(ai , f (aj )),∀i, ∀jとなる。これを成分レベルで表現すると

hAei , ej i = hei , Aej i,∀i, ∀jhAei , ej i = hei , t Aej i,∀i, ∀jとなる。一方、常にであり、上と比較することにより、

A = tA

行列計算における数式処理ソフト maxima の利用について

89を得る。すなわち、対称変換 f の任意の正規直交基底に関する表現行列は対称行列である。λをこの対称行列 A の固有方程式の解の一つとしよう。λ は一般には複素数値と考えられるのであるが、f が対称変換、すなわち A が対称行列ならば固有方程式の解は必ず実数になる。
なぜならば、v を Cn における A の λ に属する固有ベクトルとすると、

λ|v|2 = hλv, vi = hAv, vi = hv, Avi = λhv, vi = λ|v|2となって、λ = λ を得るからである。また、λ 6= µ を共に f の固有値とし、v, w をそれ
ぞれ λ, µ に属する固有ベクトルとすると、

λg(v, w) = g(λv, w) = g(f (v), w) = g(v, f (w)) = g(v, µw) = µg(v, w)
より、g(v, w) = 0 を得る。すなわち、f の異なる固有値に属する固有ベクトルは互いに直交する。
さて、この小節の目標は対称変換 (行列）の対角化可能性を示すことであるが、上の下線部で述べたことの他に、もう一つ必要な事柄がある。それは少し難しいので、詳細については線形代数学の本に譲ることにして、ここでは、事実として認めることにしよう。
f の固有値 λ の重複度が k であるとき、λ に属する固有空間の次元も k である。(dim Vλ =dimhλi = k.)
以上の準備の下に、具体的な f の、すなわち対称行列 A の対角化の手順を考えよう。
(STEP1) まず、A の固有方程式 |A − λIn | = 0 を解き、相異なる固有値 λ1 , . . . , λs とそれら各々の重複度 m1 , . . . , ms を求める。ただしここで、m1 + · · · + ms = n である。
(STEP2) 各固有値 λi , i ∈ {1, . . . , s} に対して、λi に属する A の固有空間 hλi i の正規直交基底を必要ならば Gram-Schmidt の正規直交化法などを使って作り、それを(1)(m )vi , . . . , vi i とする。
(STEP3) STEP2 で求めた固有ベクトルを
(1)

(m )

(1)

(m )

P = (v1 , . . . , v1 1 , v2 , . . . , v2 2 , . . . , vs(1) , . . . , vs(ms ) )(j)の順に並べる。すると、上述より P は Rn の正規直交基底である。各 vi , i ∈ {1, . . . , s}, j ∈{1, . . . , mi } は A の λi に属する固有ベクトルであるから、最初に議論したように、



0
λ1




..


.




λ


1


λ2




.
..


−1

P AP = 


λ2




..
.




λs




...








0
λs
となり、A は P によって対角化されることが結論される。また、P は正規直交基底であるから、行列としては直交行列である。

90

足立健朗

以上をまとめよう。任意の対称行列 A は、ある直交行列 P により対角化される。出てくる対角行列の各対角成分は A の固有値である。また直交行列 P は A の固有ベクトルからなる正規直交基底を行列と見なしたものである。
注意 7. 対称変換 f を対角化するということは、上の直交行列 P を求めることである。
本質的なのは A ではなく、P と対角行列 diag(λ1 , . . . , λs ) の方である。
さて、次に係数体を C にして考える。C 上の正定値内積空間 (V, g) と、V 上の線形変換 f : V → V についても R 上での場合とほぼ同様の議論ができる。f がエルミート変換
であるとは、条件
g(f (v), w) = g(v, f (w)), ∀v, ∀w ∈ Vが成り立つことである。このとき、V の任意の正規直交基底 A に対して、f の基底対
(A, A) に関する表現行列が、条件
A = tAを満たすことは容易にチェックできる。すなわち、エルミート変換の任意の正規直交基底に関する表現行列はエルミート行列である。この場合の、上述に対応する主張とは次の通りである。エルミート変換 f (エルミート行列 A) は、あるユニタリ行列 P によって対角化される。すなわち、



0
λ1






λ1
−1

P AP = 
...








0
λn
が成り立つ。この場合も、固有値 λ1 , . . . , λn はすべて実数である。詳細については、最早
ここでは述べない。線形代数学の本を参照して欲しい。以下の例としては取り上げることにしよう。
例 73. 実対称行列



− 21
22

2
A =  11
7
22

2
11
3
− 11
14
11



7
22
14 
11
27
22

を上の手順に従って対角化しよう。
(C1) load("gson.mac");
Warning - you are redefining the MACSYMA function EIGENVALUES
Warning - you are redefining the MACSYMA function EIGENVECTORS(D1) gson.mac(C2) 
A : matrix([-21/22,
2/11, 7/22], [2/11, -3/11, 14/11],[7/22, 14/11, 27/22]);

2
7
− 21
22
11
22
2
3
14 
− 11
(D2)  11
11
7
22

14
11

27
22

(C3) 
I : ident(3);
1 0 0(D3)  0 1 0 
0 0 1

行列計算における数式処理ソフト maxima の利用について

91(C4) ratsimp(determinant(A - x * I));
(D4) − x3 + 3x + 2(C5) chply : ratsimp(expand(charpoly(A, x)));(D5) − x3 + 3x + 2(C6) solve([chply], [x]);(D6) [x = 2, x = −1](C7) eigenvalues(A);(D7) [[2, −1] , [1, 2]](C8) 
xx : transpose(matrix([s, t, u]));s(D8)  t u(C9) 
eq1 : (A - 2*I) 
. xx;
7u
2t
65s
+ 11 − 22
22
2s 
− 25t
+ 11
(D9)  14u
11
11
7s
− 17u
+ 14t
+ 22
22
11
(C10) solve([eq1[1][1], eq1[2][1], eq1[3][1]], [s, t, u]);Dependent
eliminated: ¤¤ (3)££ equations
4%R1
,
t
=
, u = %R1
(D10)
s = %R1
77(C11) eq2 : (A - (-1)*I)
. xx;

7u
2t
s
+
+
22
11
22
8t
2s 
+
+
(D11)  14u
11
11
11
49u
14t
7s
+ 11 + 22
22(C12) solve([eq2[1][1], eq2[2][1], eq2[3][1]], [s, t, u]);Dependent equations eliminated: (2 3)(D12) [[s = −4%R3 − 7%R2, t = %R3, u = %R2]](C13) £eigenvectors(A);£¤ £¤¤(D13) [[2, −1] , [1, 2]] , [1, 4, 7] , 1, 0, − 17 , 0, 1, − 74(C14) PP : eigmat(A);

1 1
0
1 
(D14)  4 0
1
7 − 7 − 47(C15) P : gson(PP);
 √1

7
√
√
− 5252
66
5 2
33
√5

0
(D15)  √466
33
7
7
527
√
− 57√2 − 52 √33
66(C16) D : ratsimp(P^^(-1). A . P);
2 00
(D16)  0 −1 0 
0 0 −1
この例では、固有値 −1 は重解であり、それに属する固有空間は 2 次元である。行列 A
は直交行列 P 、すなわち単位ベクトルからなる固有ベクトルたちを並べたもので対角化
される。
¤

92

足立健朗

上のようにすると、途中の経過がわかるが、手っ取り早く結果のみ欲しい場合のために、専用のコマンドを作成しよう。次のプログラムを “gson.mac” に追加する。

dzmat(A) := ratsimp(gson(eigmat(A)))$
diagonalize(A) := block([B, C],B : dzmat(A),C : B^^(-1) . A . B,
ratsimp(C))$
ここで、定義したコマンド diagonalize は引数の行列を対角化した行列を返す。すなわち、引数の行列の固有値を対角成分に持つ対角行列を値として返す。また、コマンド dzmat
は引数の行列を対角化するのに使う直交行列を値として返す。
例 74.
(C1) load("gson.mac");
Warning - you are redefining the MACSYMA function EIGENVALUES
Warning - you are redefining the MACSYMA function EIGENVECTORS(D1) gson.mac(C2) 
A : matrix([2,-1, 1], [-1, 3, 2], [1, 2, -1]);
2 −1 1

−1 3
2 
(D2)
1
2 −1
(C3) 
P : dzmat(A);

√1
√ √2
√ 2 √
11
9 5+25
25−9
5
√
√


5−3
3 
 √ 5+3√(D3)  2 9√5+25 − 2 25−9√5 − √11 √√
√3 5−5√− √3 √5+5− √1
2

9 5+252(C4) 
diagonalize(A);− √55 0 0√5
0 (D4)  025−9 5

11

5

0
0 4
(C15) ratsimp(P^^(-1)
 . A . P);− √55 0 0√5
0 (D5)  05
0
0 4
エルミート行列の例についても見ておこう。
例 75.
(C1) load("gson.mac");
Warning - you are redefining the MACSYMA function EIGENVALUES
Warning - you are redefining the MACSYMA function EIGENVECTORS(D1) gson.mac(C2) A : matrix([2, 1+3*%i], [1-3*%i, -1]);

行列計算における数式処理ソフト maxima の利用についてµ¶
23i + 1(D2)
1 − 3i −1(C3) chply : ratsimp(charpoly(A , x));(D3) x2 − x − 12(C4) solve([chply],[x]);(D4) [x = −3, x = 4]
(C5) µI : ident(2);
¶
1 0
(D5)
0 1
(C6) µxx :¶transpose(matrix([s, t]));s(D6)t(C7) µeq1 : (A - (-3)*I). xx;¶(3i + 1) t + 5s
(D7)
2t + (1 − 3i) s(C8) solve([eq1[1][1], eq1[2][1]], [s, t]);Dependent
equations eliminated:
(2)
££
¤¤
(D8)
s = − 3i%R1+%R1
,
t
=
%R1
5
(C9) µeq2 : (A - 4*I)¶ . xx;(3i + 1) t − 2s
(D9)
(1 − 3i) s − 5t(C10) solve([eq2[1][1], eq2[2][1]], [s, t]);Dependent
eliminated:(2)¤¤££ equations
,
t
=
%R2
(D10)
s = 3i%R2+%R22(C11) eigenvalues(A);(D11) [[−3, 4] , [1, 1]](C12) £eigenvectors(A);£¤ £¤¤, 1, − 3i−1
(D12) [[−3, 4] , [1, 1]] , 1, 3i−1
25(C13) µeigmat(A); ¶
1
1
(D13)
3i−1
3i−1
− 5
2
(C14) Ãdzmat(A);
!
√
√
(D14)

√2
√ 7√
3 2i−
√ 2
2 7

√5
√ 7 √
√ 5
− 3 55i−
7(C15) µdiagonalize(A);
¶
−3 0
(D15)
0 4

93

94

足立健朗

10.3. 行列の巾乗再論. 行列の対角化を利用して、対角化可能な行列の巾乗を簡単に計算することができる。正方行列 A が、正則行列 P によって、


0
a1






a2
−1


P AP = 
...







0
an
の形に対角化されたとしよう。このとき、(P −1 AP )k = (P −1 AP )(P −1 AP ) · · · (P −1 AP )= P −1 Ak P




k
0
a1




k


a2

=
...








0
ak

を得る。したがって、

Ak = P (P −1 Ak P )P −1




k
0
a1




k

 −1
a2
P
=P
.


..






0
ak

を得る。これは、定義通りに巾乗を計算するのに比べて、k が大きいときは計算量を大幅に節約することができる計算法である。また、これにより任意の k に対する公式を与えていることも注目に値する。前小節で扱ったように、実対称行列あるいはエルミート行列
は直交行列あるいはユニタリ行列で対角化可能であったから、これらの行列に対してはこの方法で巾乗を計算することができる。
例 76.
(C1) load("gson.mac");
Warning - you are redefining the MACSYMA function EIGENVALUES
Warning - you are redefining the MACSYMA function EIGENVECTORS(D1) gson.mac(C2) A : matrix([2, -1, 1], [-1, 3, 2], [1, 2, -1]);

行列計算における数式処理ソフト maxima の利用について
2 −1 1
2 (D2)  −1 3
1
2 −1
(C3) 
P : dzmat(A);

√1
√ √2
√ 2 √
11
9 5+25
25−9
5
√
√


5−3
3 
 √ 5+3√(D3)  2 9√5+25 − 2 25−9√5 − √11 √√
√3 5−5√− √3 √5+5− √1112

9 5+252

25−9 5(C4) 
B : diagonalize(A);− √55 0 0√5
0 (D4)  05

0
0 4(C5) 
B5 : B^^5;

√
− 125
0
0
5
125
√
0 
(D5)  05
0
0 1024(C7) A5 : ratsimp(P . B5 . P^^(-1))(C8) 
A5;
134 −277 −59(D8)  −277 831 302 −59 302
59

95

96

足立健朗

11. ベクトル解析
11.1. 外積. この節は付録的な節である。本格的にベクトル解析を maxima で扱おうと言う意図は全くない。ベクトル解析の初歩のところで出てくる空間ベクトルの内積, 外積と、３種類の一階微分作用素 grad、div、rot に限って、これらを maxima で電卓的に扱う
方法 (の一つを) 解説する。これらは、本格的な使用には耐えないかもしれないが、ちょっとした計算には活用できると期待する。
この節では、空間ベクトルを 3 次の実列ベクトルとして扱う。まず、この小節では２つの空間ベクトルの内積と外積を、２項演算子として扱う。f が、ある集合 U 上の２項演算であるというのは、f が U × U から 値の集合 V への写像であることに他ならない。しかし、(中置きの）演算という言葉を使う場合は、値を通常の用に f (a, b) のように書くのではなくて、af b のように表す。例えば、加法 “+” なども２項演算の例である。
空間ベクトルの内積は、すでに扱った 3 次元列ベクトル空間 R3 の標準内積と同じものであるが、ここでは関数 (maxima では コマンド inprod で計算できた) ではなく、２項演算として実現しよう。maxima で新たに中置きの２項演算子を定義するには、コマン
ド infix を使う。新たに f という名の演算子を定義したい場合には、“infix(f )” で f が新しい演算子であることを宣言し、その後 “a f b := 関数定義式” で実体を定義する。ここで、a と b は仮引数である。それでは、空間ベクトルの内積を @∗ という記号で表すことにして、maxima で使うために定義しよう。次のプログラムを “vectanal.mac” と言うファイルに保存する。
infix("@*")$
v @* w := block([val],
val : v[1][1] * w[1][1] + v[2][1] * w[2][1] + v[3][1] * w[3][1] ,return(val))$
早速、実験してみよう。
例 77.
(C1) load("vectanal.mac");(D1) vectanal.mac(C2) 
v : matrix([a],[b],[c]);a(D2)  b c(C3) 
w : matrix([s],[t],[u]);s(D3)  t u(C4) v @* w;(D4) cu + bt + as(C5) v1 : transpose(matrix([1, -2, 3]));

行列計算における数式処理ソフト maxima の利用について

97




1
(D5)  −2 
3(C6) 
w1 : transpose(matrix([4,1, -2]));
4(D6)  1 −2(C7) v1 @* w1;(D7) − 4
内積というのは２つの空間ベクトルに一つの実数を対応させる２項演算子であるが、３次元、したがって空間ベクトルの場合には、外積と言うもう一つの有用な２項演算子がある。外積は２つの空間ベクトルに一つの空間ベクトルを対応させるもので、その定義は座標を使って次のように定義される。
    
a1b1
a2 b3 − b2 a3
a2  ×  b2  = a3 b1 − a1 b3  .
a3
b3
a1 b2 − a2 b1
ここで、“×” が数学で普通外積を表すのに使われる記号である。この記号は、キーボー
ド上に存在しないので maxima で外積を表す記号としては、“@x” を使うことにしよう。
内積のときのように定義する。以下のプログラムをファイル “vectanal.mac” に追加する。

infix("@x")$
v @x w := block([valvec],
valvec : zeromatrix(3,1),
valvec[1][1] : v[2][1] * w[3][1] - v[3][1] * w[2][1] ,valvec[2][1] : v[3][1] * w[1][1] - v[1][1] * w[3][1] ,valvec[3][1] : v[1][1] * w[2][1] - v[2][1] * w[1][1] ,return(valvec))$
これで、maxima で外積を使う準備ができた。
例 78.
(C1) load("vectanal.mac");(D1) vectanal.mac(C2) 
v : transpose(matrix([a,b, c]));
a(D2)  b c(C3) 
w : transpose(matrix([s,t, u]));
s(D3)  t u(C4) v @x w;

98

足立健朗




bu − ct(D4)  cs − au at − bs(C5) 
v1 : transpose(matrix([1,-2, 3]));
1(D5)  −2 
3(C6) 
w1 : transpose(matrix([4,1, -2]));
4(D6)  1 −2(C7) 
v1 @xw1;
1

14 
(D7)
9
外積の幾何学的意味は次の通りである。v、w を共にゼロベクトルではなく、互いに平行ではない２つの空間ベクトルとする。このとき v と w の外積 v × w は、v と w の両方に対して垂直な空間ベクトルで、その長さは v と w の張る平行四辺形の面積 |v||w|| sin θ|に等しい。ここで、θ は v と w の間の角とする。これらのことを maxima でチェックしてみよう。確かめればいいのは、v • (v × w) = w • (v × w) = 0 と、

|v × w|2 = |v|2 |w|2 sin2 θ= |v|2 |w|2 (1 − cos2 θ)22
2
2 |v| |w| − (v • w)
= |v| |w|2
2
2

|v| |w|2

2= |v| |w| − (v • w)2
である。
例 79.
(C1) load("vectanal.mac");(D1) vectanal.mac(C2) 
v : transpose(matrix([a,b, c]));
a(D2)  b c(C3) 
w : transpose(matrix([s,t, u]));
s(D3)  t u(C4) v @* (v @x w);(D4) a (bu − ct) + b (cs − au) + c (at − bs)(C5) ratsimp(D4);(D5) 0(C6) w @* (v @x w);

行列計算における数式処理ソフト maxima の利用について

99(D6) s (bu − ct) + t (cs − au) + (at − bs) u(C7) ratsimp(D6);(D7) 0(C8) ratsimp((v @x w) @* (v @x w) - ((v @* v) * (w @* w) - (v @* w)^2));(D8) 0外積というのは、非可換な積である。一般に v × w = −w × v が成り立つ、すなわち、積の順序を入れ換えると計算結果が逆向きのベクトルになる。したがって、特に v × v = 0となる。外積と内積を使って計算する量 (u × v) • w は、u、v、w がこの順序で張る平行六面体の符合つき体積に等しく、それは行列 A = (u, v, w) の行列式 |A| に等しい。これらのこともチェックしてみよう。
例 80.
(C1) load("vectanal.mac");(D1) vectanal.mac(C2) 
v : transpose(matrix([a,b, c]));
a(D2)  b c(C3) 
w : transpose(matrix([e,f, g]));
e(D3)  f g(C4) 
v @x w; 
bg − cf(D4)  ce − ag 
af − be(C5) 
- (w @x v);
bg − cf(D5)  ce − ag 
af − be(C6) 
ratsimp((v@x w) + (w @x v));0(D6)  0 
0(C7) 
v @xv;
0

0 
(D7)
0
(C8) 
w @xw;
0(D8)  0 
0(C9) kill(v, w);(D9) DONE(C10) u : transpose(matrix([a, b, c]));

100

足立健朗




a
(D10)  b c(C11) v : transpose(matrix([d,e, f]));
d(D11)  e f(C12) w : transpose(matrix([g,h, i]));
g(D12)  h i(C13) vol1 : ratsimp((u @x v) @* w);(D13) (ae − bd) i + (cd − af ) h + (bf − ce) g(C14) B : zeromatrix(3,3)(C15) B[1] : transpose(u)[1];(D15) [a, b, c](C16) B[2] : transpose(v)[1];(D16) [d, e, f ](C17) B[3] : transpose(w)[1];(D17) [g, h, i](C18) A : transpose(B);a d g(D18)  b e h 
c f i(C19) vol2 : ratsimp(determinant(A));(D19) (ae − bd) i + (cd − af ) h + (bf − ce) g(C20) ratsimp(vol1 - vol2);(D20) 011.2. 微分作用素. ３次元空間 R3 上のベクトル場というのは、空間の各点にひとつのベクトルを対応させる規則のことである。R3 から、あるベクトル空間 V への写像のことと言ってもいい。特に、R3 の各点 x に x を起点とする有向線分の一つを対応させるものは R3 上の接ベクトル場と呼ばれるが、これは一つの写像 f : R3 → R3 と思っていい。
この小節では、３次元空間 R3 に正規直交系として標準基底 E = (e1 , e2 , e3 ) を入れて考え、空間内の点を E に関する座標 x = t (x, y, z) をつかって表す。通常の関数と区別する
ため、ベクトル場は太字 f などを使って表すことにしよう。f を接ベクトル場とする。fの点 x = t (x, y, z) に置ける値は R3 の一つのベクトルであり、f (x, y, z)
f (x) =  g(x, y, z) h(x, y, z)の形に書ける。ここで、f, g, h は各座標関数 x, y, z を変数とする３変数の関数である。
f が (x, y, z) に関して連続であるとき、f は連続であると言う。また、各成分 f, g, h が行列計算における数式処理ソフト maxima の利用について

101

x, y, z の各々について偏微分可能であるとき、f はその座標関数について偏微分可能であると言う。以下では、各成分が座標関数について適当な微分可能性を持っているようなベクトル場を扱う。
一階線形微分作用素 grad は、R3 上の関数 f に作用して、R3 上の接ベクトル場 grad fを返す作用素である。ここで、f は各座標関数について偏微分可能な関数で、grad f は
 ∂f

(x)
∂x
(grad f )(x) =  ∂f (x)
∂y
∂f
(x)
∂z

により定義される接ベクトル場である。逆に、一階線形微分作用素 div は、各座標関数
について偏微分可能な接ベクトル場


f (x)
f (x) =  g(x) h(x)に対して、関数 div を対応させる作用素である。ここで、 div は(div f )(x) =

∂f
∂g
∂h
(x) +
(x) +
(x)
∂x
∂y
∂z

で定義される。もう一つの一階線形微分作用素 rot は、各座標関数について偏微分可能
な接ベクトル場
f (x)f (x) =  g(x) h(x)に対して、接ベクトル場 rot f を対応させる作用素である。ここで、rot f は
´
³

∂g
∂h
−
(x)
∂y
∂z
 ¡ ∂f¢∂h(rot f )(x) = ³ ∂z − ∂x ´ (x) 

∂g
∂f
−
(x)
∂x
∂y
で定義される。grad f は関数 f の勾配ベクトル場、div f はベクトル場 f の発散、rot f はベクトル場 f の回転という。
これらの微分作用素は、古典電磁気学の完成形である Maxwell の理論の記述に使われることからもわかるように、極めて基本的で重要な微分作用素と考えられる。詳しい
数学的性質や、その物理現象への応用などは相応の専門書に譲ることにして、ここでは、
maxima でこれらの作用素を、前置演算子として定義し、いくつかの計算例を示すことで満足することにしよう。
maxima で前置演算子を定義するのには、コマンド prefix を使う。prefix の詳細はマニュアルを見てもらうことにして、大体のところは infix と同様につかえると思えば良い。
maxima で関数の微分を扱うコマンドは diff である。diff(f, x); とタイプすると、関数 fを x に関し (偏) 微分したものを返す。ただし、f は具体的に x を含む形で定義された関数か、あるいは depends(f, x) というコマンドで、x に依存することを宣言された形式的な関数である必要がある。それ以外のときは、x に関しては定数とみなされ、0 を返す。
prefix と diff を使って、コマンド grad、div および rot を定義しよう。次のプログラム
をファイル “vectanal.mac” に追加する。

102

足立健朗

prefix("Grad", 142, expr, expr)$
prefix("Div", 142, expr, expr)$
prefix("Rot", 142 , expr, expr)$"Grad"(f) := block([gv],gv : zeromatrix(3,1),
gv[1][1] : diff(f, x),
gv[2][1] : diff(f, y),gv[3][1] : diff(f, z),return(gv))$"Div"(v) := block([f],f : diff(v[1][1], x) + diff(v[2][1], y) + diff(v[3][1], z),return(f))$"Rot"(v) := block([w],w : zeromatrix(3,1),
w[1][1] : diff(v[3][1], y) - diff(v[2][1], z),w[2][1] : diff(v[1][1], z) - diff(v[3][1], x),w[3][1] : diff(v[2][1], x) - diff(v[1][1], y),return(w))$
Laplacian(f) := div (grad (f))$
いくつかの使用例を見てみよう。
例 81. noindent(C1) load("vectanal.mac");(D1) vectanal.mac(C2) diff(f, x);(D2) 0(C3) depends(f , x);(D3) [f (x)](C4) diff(f, x);df(D4) dx(C5) depends(f, [x, y, z]);(D5) [f (x, y, z)](C6) 
grad f;
dfdx(D6)  d f 
dy
d
f
dz

行列計算における数式処理ソフト maxima の利用について

103(C7) 
F : transpose(matrix([o, p, q]));o(D7)  p q(C8) div F;(D8) 0(C9) 
rot F;0(D9)  0 
0(C10) depends([o, p, q], [x, y, z]);(D10) [o (x, y, z) , p (x, y, z) , q (x, y, z)](C11) div F;ddd(D11) dz
q + dyp + dxo(C12) rot F;dd
q − dzpdydd(D12)  dz
o − dxq ddp − dy odx(C13) div (grad f);
d2
d2
d2(D13) dz
2 f + dy 2 f + dx2 f(C14) grad (div F);

d2
d2
d2
q
+
p
+
o
dxdz
dxdy
dx2
 d2

d2
d2
q + dyp+o (D14)  dydz
2
dxdy
d2
d2
d2
q + dydz
p + dxdzodz 2(C15) rot 
(grad f);0(D15)  0 
0(C16) div (rot F);(D16) 0(C17) Laplacian(f);d2d2d2(D17) dz
2 f + dy 2 f + dx2 f
div(gradf ) は通常 ∆(f ) あるいは ∆f などと書かれ、f のラプラシアンと呼ばれる。ラプラス作用素 ∆ は２階の楕円型作用素の一つである。maxima のコマンドとしては Laplacianとして使える。
これらの微分作用素が実際にどのように使用されるかについては、ここでは全く解説しないので、適当な物理の本を読んでいただきたい。何を計算しているかについては説明しないが、最後にもうちょっとだけ、計算例を付け加えよう。
例 82.
(C1) load("vectanal.mac");(D1) vectanal.mac(C2) g(x,y,z) := -1/sqrt(x^2 + y^2 + z^2);(D2) g (x, y, z) := √ 2 −12 2
x +y +z

104

足立健朗

(C3) 
grad g(x,y,z);x3(z 2 +y 2 +x2 ) 2



y
3 
(D3) 
 (z2 +y2 +x2 ) 2 z

3(z 2 +y 2 +x2 ) 2(C4) ratsimp(sqrt(D3 @* D3));1(D4) √ 4
2
2 2
4
2 2
4
z +(2y +2x )z +y +2x y +x(C5) factor(z^4 + (2*y^2 + 2*x^2)*z^2 + y^4 + 2*x^2*y^2 + x^4);2(D5) (z 2 + y 2 + x2 )(C6) ratsimp(D4^2 - 1 / D5);
(D6) 0(C7) pr : sqrt(x^2 + y^2 + z^2);(D7)
z 2 + y 2 + x2
(C9) 
rot (grad g(x, y, z));0(D9)  0 
0(C10) div (grad g(x,y,z));222(D10) 2 23 2 3 − 2 3z2 2 5 − 2 3y2 2 5 − 2 3x2 2 5(z +y +x ) 2(z +y +x ) 2(z +y +x ) 2(z +y +x ) 2(C11) ratsimp(D10);(D11) 0(C12) F : transpose(matrix([sin(x*y),cos(x*y), z^2 - x*y]));sin (xy)(D12)  cos (xy) z 2 − xy(C13) div F;(D13) 2z − x sin (xy) + y cos (xy)(C14) rot F;
−xy(D14) −y sin (xy) − x cos (xy)(C15) rot (rot F);x2 sin (xy) − sin (xy) − xy cos (xy)(D15)  −xy sin (xy) + y 2 cos (xy) + cos (xy) 0(C16) grad (rot F);

行列計算における数式処理ソフト maxima の利用について

 

105

 
−1
 
 
0


2


xy
sin
(xy)
−
y
cos
(xy)
−
cos
(xy)
 
 


0



 
1
(D16) 
x2 sin (xy) − sin (xy)
 − xy cos (xy)


0



 0 

0
(C17) div (rot F);(D17) 0(C18) G : 1 / r^2;(D18) z2 +y12 +x2
上の (C16) の grad(rot F ) は数学的には定義できないが、maxima での演算子の特性として rot F の各成分に対して grad を適用したものが値として返されている。

106

足立健朗

12. あとがきというか言い訳
大学の線形代数の授業の補助として、行列計算の算盤がわりとして maxima を使うための手軽な入門書を目指して書き始めた本文書ですが、当初の意に反して、マニアックな記述に満ちた、著者自身のための備忘録となってしまいました。そういうわけで、この文書を読んで失望される方もおられると思いますが、只々、お詫び申し上げるのみです。
この文書自体は、maxima の入門書でも、線形代数学の入門書でもありません。それらについては、それ用の文献を見て下さい。この下に、私が好きな参考書のリストを掲載します。
この文書を書くにあたっては、私の３人の友人に読者、モニターになって頂き、貴重な御意見を頂きました。上田秀雄さん、古閑義之さん、吉岡雄一さんに感謝申し上げます。
maxima 関係の参考資料[1] Maxima Manual, http://maxima.sourceforge.net/docs.shtml .
[2] Paulo Ney de Souza, Richard J. Fateman, Joel Moses and Cliff Yapp,
The Maxima Book, http://maxima.sourceforge.net/docs.shtml .
[3] Michael Clarkson, DOE-Maxima Reference Manual,
http://maxima.sourceforge.net/docs.shtml .
線形代数学の参考文献[4] 伊原信一郎, 河田敬義, 線型空間、アフィン幾何, 岩波基礎数学選書[5] 杉浦光夫, 横沼健雄, ジョルダン標準形、テンソル代数, 岩波基礎数学選書[6] 田坂隆士, ２次形式, 岩波基礎数学選書[7] 斎藤正彦, 線型代数入門, 東京大学出版会[8] 長岡亮介, h 新版 i 線型代数学-入門と展望, ブレーン出版
その他の参考文献[9] Joachim von zur Gathen and Jürgen Gerhard,
Modern Computer Algebra, Cambridge University Press.
[10] 酒井皇治, やさしい Lisp 入門, 株式会社ニコラステラ[11] ポール グレアム, ANSI Common Lisp, ピアソン • エデュケーション[12] Paul Graham, On Lisp, http://www.paulgraham.com[13] David S. Touretzky, Common Lisp: A Gentle Introduction to Symbolic Computation, http://www-2.cs.cmu.edu/ dst/LispBook/index.html
E-mail address: kenrou@yo.rim.or.jp