線形代数の演習問題と解答の作成

1.乱数行列の発生

指定した範囲の乱整数を使った乱数行列を発生する方法。例えば−7から9までの整数を要素とする乱数行列は、

> f:=rand(-7..9):m:=2;n:=3;

m := 2

n := 3

> A:=matrix(m,n,[seq(f(),k=1..m*n)]);

A := matrix([[-1, -7, 1], [3, 9, -7]])

または、次のようにしてm行n列の乱数行列が発生できる。

> f:=rand(-7..9):

> A:=(m,n)->matrix(m,n,[seq(f(),i=1..m*n)]);

A := proc (m, n) options operator, arrow; matrix(m,...

> A(3,4);

matrix([[7, 6, -4, -1], [2, -7, -1, -5], [6, 6, 4, ...

またfor loopを使って乱数行列を望みの数だけ発生できる。

> for k to 3 do: A(3,3): od;

matrix([[9, 3, -4], [8, 9, -5], [-7, -6, -6]])

matrix([[-3, -1, -3], [0, 1, 0], [-3, -2, 9]])

matrix([[-6, -7, -3], [-3, 3, -1], [-7, 9, 6]])

>

2.行列の演算
2.1 行列の和と差
行列の和と差は同じ形の行列同士でなければならないことを考慮して、

> restart:

> WA:=proc(m,n)

> local f,A,B,C:

> f:=rand(-7..9):

> A:=matrix(m,n,[seq(f(),i=1..m*n)]):

> B:=matrix(m,n,[seq(f(),i=1..m*n)]):

> C:=evalm(A+B):

> print(matrix(A)+matrix(B)=matrix(C));

> end:

> WA(3,4);

matrix([[-1, -7, 1, 3], [9, -7, 7, 6], [-4, -1, 2, ...

>

 なおendの前の2行を,
> C:=evalm(A-B):
> print(A,` - `,B,` = `,C);
と変更すれば行列の差の演習問題が作れる。また,連続して多数の演習問題を作りたいときは、
> for j to 3 do: WA(3,2);od;

>

2.2 行列の積
行列の積の場合は(m×n)と(n×p)行列の積のみが定義されている。

> SEKI:=proc(m,n,p)

> local f,A,B,C:

> f:=rand(-7..9):

> A:=matrix(m,n,[seq(f(),i=1..m*n)]):

> B:=matrix(n,p,[seq(f(),i=1..n*p)]):

> C:=evalm(A&*B):

> print(matrix(A) * matrix(B) = matrix(C));

> end:

たとえば,(3×2)行列と(2×5)行列の積の演習問題と答は、

> SEKI(3,2,5);

>

matrix([[-7, -6], [-6, -3], [-1, -3]])*matrix([[0, ...

>

2.3 逆行列
まずwith(linalg)で線形代数パッケージを呼び出しておく。
行列Aの逆行列はinverse(A)という組み込みコマンドで求まるが、これでは答えの逆行列の要素が分数となり面白くない。逆行列の各要素も整数となる演習問題を作成するプロシージャを示すと、

> restart:

> gyaku:=proc(n::{2,3})

> local f,A,B,d,i,j:
with(linalg):

> f:=rand(-7..9):i:=0:
for j to 2000 while i<3 do

> A:=matrix(n,n,[seq(f(),k=1..n^2)]):

> d:=det(A):
if abs(d)=1
then i:=i+1:
print( cat(`No. `,i),` A` = matrix(A),` inv A`=inverse(A),` det(A)`=d) fi:

> od:

> end:

>

このプロシージャは,与えられた行列の行列式の値が1または−1であれば逆行列の全要素も整数となることを利用する。2000回トライし求める逆行列を最大3個まで書き出す。このような行列は数百回に1度程度の割合で発生するので、2000回のトライ中3個発生しない場合もある。

>

> gyaku(3);

`No. 1`, ` A` = matrix([[-3, 4, 1], [5, -6, -2], [-...

`No. 2`, ` A` = matrix([[3, 8, 6], [1, 0, -1], [-4,...

`No. 3`, ` A` = matrix([[-5, 4, 4], [-7, 5, 7], [0,...

2.4 余因子行列
余因子行列を求めるコマンド adjoint(A) を利用し、

> with(linalg):f:=rand(-7..9):n:=3:

> A:=matrix(n,n,[seq(f(),i=1..n^2)]);

A := matrix([[-3, 8, 2], [-1, 2, 7], [3, 1, -7]])

> adjoint(A);

matrix([[-21, 58, 52], [14, 15, 19], [-7, 27, 2]])

2.5 行列のランク 
行列のランクは殆ど行列の次数と一致する。行列の次数とランクが異なる例題を見つけるのプロシージャを示す。4次の行列でランクが3以下の乱数行列は1000回に1回程度しか現れない。

> restart:with(linalg):

Warning, the protected names norm and trace have been redefined and unprotected

> ranking:=proc(n)

> local f,k,A,j,jj:

> if nargs=2 then k:=args[2]:else k:=1:fi:

> f:=rand(-7..9):jj:=0:

> for j to 3000 while jj<k do

> A:=matrix(n,n,[seq(f(),i=1..n^2)]):

> if rank(A)<n then jj:=jj+1:

print(`A`=matrix(A),` rank(A)`=rank(A));fi:

> od:

> end:

> ranking(3);

A = matrix([[2, 1, -4], [9, -6, 3], [2, -4, 6]]), `...

> ranking(4);

A = matrix([[-7, -5, -5, 6], [-7, 3, -5, 6], [4, 5,...

なお > ranking(4,3);とすれば4次の行列を3個出力するが、3000回のトライ中3個生成されないときもある。

>

3 行列式とその展開
3.1 行列式の値
n次の行列とその行列式の値を求める演習問題を作成する。determ(n,k)と第2の引数を指定すると連続してk個の演習問題が出力される.

> restart:with(linalg):

Warning, the protected names norm and trace have been redefined and unprotected

> determ:=proc(n)

> local f,k,A,j:

> if nargs=2 then k:=args[2]:else k:=1:fi:

> f:=rand(-7..9):

> for j to k do

> A:=matrix(n,n,[seq(f(),i=1..n^2)]):

> print(`A`=matrix(A),` det(A)`=det(A));

> od:

> end:

> determ(4,2);

A = matrix([[-1, -7, 1, 3], [9, -7, 7, 6], [-4, -1,...

A = matrix([[4, 6, 9, 3], [-4, 8, 9, -5], [-7, -6, ...

>

3.2 行列式の展開
任意の乱数行列を与え,それを指定した行または列で展開して行列式の計算を行う。

行による展開

> restart:with(linalg):

Warning, the protected names norm and trace have been redefined and unprotected

> gyouten:=proc(A,k)

> local i,gk,d0,d1,d2,n:

> n:=rowdim(A):
if k>n then RETURN(`ERROR, k is larger than `,n) fi:

> d0:=add((-1)^(i+k)*A[k,i]*minor(A,k,i),i=1..n);

> d1:=[seq((-1)^(i+k)*A[k,i]*det(minor(A,k,i)),i=1..n)];

> d2:=add((-1)^(i+k)*A[k,i]*det(minor(A,k,i)),i=1..n);

> print(A=d0);

> print(d1,`det(A)`=d2);

> end:

>

列による展開

> retsuten:=proc(A,k)

> local i,gk,d0,d1,d2,n:

> n:=coldim(A):
if k>n then RETURN(`ERROR, k is larger than `,n) fi:

> d0:=add((-1)^(i+k)*A[i,k]*minor(A,i,k),i=1..n);

> d1:=[seq((-1)^(i+k)*A[i,k]*det(minor(A,i,k)),i=1..n)];

> d2:=add((-1)^(i+k)*A[i,k]*det(minor(A,i,k)),i=1..n);

> print(A=d0);

> print(d1,`det(A)`=d2);

> end:

n次の乱数行列の発生

>

> f:=rand(-7..9):n:=5:

> A:=matrix(n,n,[seq(f(),k=1..n^2)]);

A := matrix([[-1, -7, 1, 3, 9], [-7, 7, 6, -4, -1],...

Aの3行目による展開

> gyouten(A,3);

A = 2*matrix([[-7, 1, 3, 9], [7, 6, -4, -1], [4, 6,...

[9478, -13041, -338, 27330, 22560], `det(A)` = 4598...

>

Aの2列目による展開

> retsuten(A,2);

A = 7*matrix([[-7, 6, -4, -1], [2, -1, -5, 6], [6, ...

[43365, 52689, -13041, 12128, -49152], `det(A)` = 4...

>

4 連立1次方程式の解
4.1 ベクトル解
整数解をもつ連立方程式を作成する。

> restart:

> sim:=proc(n)

> local A,X,b,i,j,k,x,f,s,eq:
if nargs=2 then k:=args[2]:else k:=1: fi:

> with(linalg):

> f:=rand(-7..9):
for i to k do
print(` `):

> A:=matrix(n,n,[seq(f(),k=1..n^2)]):

> X:=matrix(n,1,[seq(f(),k=1..n)]):

> b:=multiply(A,X):

> for j to n do

> if n=2 then eq:=A[j,1]*x[1]+A[j,2]*x[2]=b[j,1]:

> elif n=3 then eq:=A[j,1]*x[1]+A[j,2]*x[2]+A[j,3]*x[3]=b[j,1]:
else eq:=A[j,1]*x[1]+A[j,2]*x[2]+A[j,3]*x[3]+A[j,4]*x[4]=b[j,1]:

> fi:

> print(eq):

> od:

> print(cat(`No.`,i),` sol`=evalm(X)):
od:

> end:

4元1次連立方程式の問題と解答は、

> sim(4);

Warning, the protected names norm and trace have been redefined and unprotected

` `

-x[1]-7*x[2]+x[3]+3*x[4] = -28

9*x[1]-7*x[2]+7*x[3]+6*x[4] = 75

-4*x[1]-x[2]+2*x[3]-7*x[4] = -25

-x[1]-5*x[2]+6*x[3]+6*x[4] = 38

`No.1`, ` sol` = matrix([[4], [6], [9], [3]])

>

4.2 クラーメルの公式

解きたい連立方程式を係数行列と定数ベクトルで表した行列で与える。クラーメルの公式で解き、分母、分子の値と方程式の解を求める。

> restart:with(linalg):

Warning, the protected names norm and trace have been redefined and unprotected

> A:=matrix([[8,1,1,7,13],[-7,0,-7,3,84],[-4,7,3,-1,57],[5,1,-1,-2,-27]]);

A := matrix([[8, 1, 1, 7, 13], [-7, 0, -7, 3, 84], ...

> n:=rowdim(A):Ad0:=det(submatrix(A,1..n,1..n)):

> for j to n do

> Ad:=det(submatrix(swapcol(A,j,n+1),1..n,1..n)):ans:=Ad/Ad0:
#printf(" %a%s %d%s%d%s %d\n",X[j],"=",Ad,"/",Ad0,"=",ans);

print(cat(`X[`,j,`] = `,Ad,` / `,Ad0,` = `,ans));

> od:

`X[1] = 17510 / -3502 = -5`

`X[2] = -28016 / -3502 = 8`

`X[3] = 14008 / -3502 = -4`

`X[4] = -24514 / -3502 = 7`

>

4.3 掃き出し法
 連立方程式における係数行列と定数列ベクトルを並べた拡大係数行列をもとに,掃き出し法で連立方程式の解を求める。

> restart:with(linalg):

Warning, the protected names norm and trace have been redefined and unprotected

> n:=4;

n := 4

> A:=matrix(n,n+1,[8,1,1,7,13,-7,0,-7,3,84,-4,7,3,-1,57,5,1,-1,-2,-27]);

A := matrix([[8, 1, 1, 7, 13], [-7, 0, -7, 3, 84], ...

> gaussjord(A);

matrix([[1, 0, 0, 0, -5], [0, 1, 0, 0, 8], [0, 0, 1...

第4列に解が現れる.

>

4.4 掃き出し法による逆行列の計算
 行列A と同じ次数の単位行列E を並べた行列[A,E] を掃き出し法により [E,B]とすることができれば Bが Aの逆行列となる.これを利用して逆行列を求める。

> restart:with(linalg):

Warning, the protected names norm and trace have been redefined and unprotected

> n:=3:f:=rand(-7..9):

> A:=matrix(n,n,[seq(f(),k=1..n^2)]);

A := matrix([[-1, -7, 1], [3, 9, -7], [7, 6, -4]])

> E:=band([1],n);

E := matrix([[1, 0, 0], [0, 1, 0], [0, 0, 1]])

> AE:=augment(A,E);

AE := matrix([[-1, -7, 1, 1, 0, 0], [3, 9, -7, 0, 1...

> EB:=gaussjord(AE);

EB := matrix([[1, 0, 0, 3/104, -11/104, 5/26], [0, ...

右半分に逆行列が生成される.

5. 固有値と固有ベクトル
行列の工学問題への応用として重要である固有値と固有ベクトルを求める。固有値が整数となる演習問題をトライアンドエラー法で作成する。3次の正方行列では1000回に1回程度の割合で望みの行列が求る。

> restart:

> eigen:=proc(n::{2,3})

> local f,A,B,i,j,nk,EV:
with(linalg):
if nargs=2 then nk:=args[2] else nk:=1 fi:

> f:=rand(1..8):i:=0:
for j to 1000 while i<nk do

> A:=matrix(n,n,[seq(f(),k=1..n^2)]):

> B:=[eigenvalues(A)]:
if is(B[1],integer)=true and is(B[2],integer)=true
then i:=i+1:
print(cat(`No.`,i)=matrix(A),`A`=matrix(A),`eigenvalues`=B);
EV:=eigenvectors(A):
print(`eigenvectors`= EV);
fi:

> od:

> end:

2次の演習問題

> eigen(2);

`No.1` = matrix([[6, 4], [2, 8]]), A = matrix([[6, ...

eigenvectors = ([4, 1, {vector([-2, 1])}], [10, 1, ...

3次の演習問題

> eigen(3);

`No.1` = matrix([[7, 7, 7], [7, 7, 8], [2, 2, 7]]),...

eigenvectors = ([0, 1, {vector([-1, 1, 0])}], [4, 1...

なお,出力中,eigenvectors(固有ベクトル)はMapleコマンドで表示される形式である。
また,eigen(3,2) とす第2引数を指定すれば3次の正方行列が2個出力される.

> eigen(3,2);

`No.1` = matrix([[4, 4, 1], [6, 2, 1], [3, 7, 5]]),...

eigenvectors = ([3, 1, {vector([1, 1, -5])}], [-2, ...

`No.2` = matrix([[1, 4, 5], [3, 2, 5], [7, 1, 8]]),...

eigenvectors = ([-2, 1, {vector([7, 1, -5])}], [13,...

>