C8MW.html

EigenValues.mw

Matrix Inverse

Examples

まず，Mを造る．
M := <<1,2,0>|<2,3,2>|<0,2,1>>;

対応する3元連立方程式

> v := <5,4,2>;
xx:=<x,y,z>;
M.xx=v;

逆行列(Matirx Inverse)を求めて，vに作用すると，{x,y,z}が求まる．

> with(LinearAlgebra):
MatrixInverse(M).M.xx=MatrixInverse(M).v;

ガウス消去法を試せるJava

> with(Student[LinearAlgebra]):
GaussianEliminationTutor( M );
GaussianEliminationTutor( M, v );

LU decomposition

> #1
restart;
with(LinearAlgebra):

行列Mと数値ベクトルbを定義．

> #2
AA[1]:=Matrix(1..4,1..4,[[4,3,2,1],[2,5,-3,-2],
[1,-4,8,-1],[-3,2,-4,5]]);
b:=Vector([20,-5,13,9]);

行列の[1,1]対角成分からAA[2,1],[3,1],[4,1]を消去する変換行列T1を求める．

> #3
T[1]:=Matrix(1..4,1..4,[[1,0,0,0],[-1/2,1,0,0],
[-1/4,0,1,0],[-(-3/4),0,0,1]]);

T1をAA[1]に作用させる．

> #4
AA[2]:=MatrixMatrixMultiply(T[1],AA[1]);

> #5
T[2]:=Matrix(1..4,1..4,[[1,0,0,0],[0,1,0,0],
[0,-AA[2][3,2]/AA[2][2,2],1,0],[0,-AA[2][4,2]/AA[2][2,2],0,1]]);

> #6
AA[3]:=MatrixMatrixMultiply(T[2],AA[2]);

> #7
T[3]:=Matrix(1..4,1..4,[[1,0,0,0],[0,1,0,0],
[0,0,1,0],[0,0,-AA[3][4,3]/AA[3][3,3],1]]);

> #8
AA[4]:=MatrixMatrixMultiply(T[3],AA[3]);

行列を三角行列に分解するために用意されているLUDecompositionコマンド

> #9
(P,L,U):=LUDecomposition(AA[1]);

変換行列の積

> #10
T[3].T[2].T[1];

下三角行列(L)の逆行列

> #11
MatrixInverse(L);

LとUの積はもとの行列AA[0]に一致する．

> #12
L.U;

変換行列を一度に数値ベクトルbに作用させると，ガウス消去法に対応した行列になる．

> #13
T[3].T[2].T[1].b;
LUDecomposition(<AA[1]|b>,output='U');

'R'を指定すると対角化した行列，あるいは連立方程式の解

> #14
LUDecomposition(<AA[1]|b>,output='R');

Partial Pivoting

> #1
restart;with(LinearAlgebra):

> #2
AA[1]:=Matrix(1..4,1..4,[[2,-4,3,-1],[1,-2,2,1],
[1,-5,4,-3],[3,2,-2,-2]]);
b:=Vector([-2,1,-8,1]);

1列目の消去をしたところ，同時にT[1].AA[1]の2行2列目の要素が0になってしまった．

> #3
t:=-1/AA[1][1,1];
T[1]:=Matrix(1..4,1..4,[[1,0,0,0],[t*AA[1][2,1],1,0,0],
[t*AA[1][3,1],0,1,0],[t*AA[1][4,1],0,0,1]]);
T[1].AA[1];

Partial pivoting, P[1]行列を左から作用すると，2行目と3行目が入れ替わる．この操作は変数の並びは変わらず単に方程式の順番を変更しただけである．

> #4
P[1]:=Matrix(1..4,1..4,[[1,0,0,0],[0,0,1,0],
[0,1,0,0],[0,0,0,1]]);
AA[2]:=P[1].T[1].AA[1];

> #5
n:=2;
t:=-1/AA[n][n,n]:
T[n]:=Matrix(1..4,1..4,[[1,0,0,0],[0,1,0,0],
[0,t*AA[n][3,n],1,0],[0,t*AA[n][4,n],0,1]]):
AA[n+1]:=T[n].AA[n];

> #6
n:=3;
t:=-1/AA[n][n,n];
T[n]:=Matrix(4,4):
for i from 1 to 4 do
T[n][i,i]:=1;
end do:
T[n][n+1,n]:=t*AA[n][n+1,n];
AA[n+1]:=T[n].AA[n];

上三角行列にした後も，対角の規格化，ならびに上三角の非対角成分を後ろ側から消去していく操作(後退代入)を引き続き作用させることで行列Aが対角化される

> #7
S[0]:=Matrix(4,4):
for i from 1 to 4 do
S[0][i,i]:=1/AA[4][i,i];
end do:
AA[5]:=S[0].AA[4];

> #8
n:=3;
S[n]:=Matrix(4,4):
S[n][4,4]:=1;
for m from n to 1 by -1 do
S[n][m,m]:=1;
S[n][m,4]:=-AA[5][m,4];
end do:
AA[6]:=S[n].AA[5];

> #9
n:=2;
S[n]:=Matrix(4,4):
S[n][4,4]:=1:
S[n][3,3]:=1:
for m from n to 1 by -1 do
S[n][m,m]:=1;
S[n][m,3]:=-AA[6][m,3];
end do:
AA[7]:=S[2].AA[6];

> #10
n:=1;
S[n]:=Matrix(4,4):
S[n][4,4]:=1:
S[n][3,3]:=1:
S[n][2,2]:=1:
for m from n to 1 by -1 do
S[n][m,m]:=1;
S[n][m,2]:=-AA[7][m,2];
end do:
S[1].AA[7];

逆行列(Z)

> #11
Z:=S[1].S[2].S[3].S[0].T[3].T[2].T[1].P[1];
Z.AA[1];

方程式の根

> #12
Z.b;

MapleのLUDecompositionコマンドをこの行列に適用すると，置換行列(permutation matrix)Pは単位行列ではなくなる．

> (P,L,U):=LUDecomposition(AA[1]);

P.A=L.U

> P.AA[1];
L.U;

Jacobi method for solving simultaneous linear equations

> #1
restart;
n:=4:
AA:=Matrix(1..4,1..4,[[4,3,2,1],[2,5,-3,-2],
[1,-4,8,-1],[-3,2,-4,5]]);
b:=Vector([20,-5,13,9]);

#2 Jacobi method

> x0:=[0,0,0,0]:
x1:=[0,0,0,0]:
for iter from 1 to 20 do
for i from 1 to n do
x1[i]:=b[i];
for j from 1 to n do
x1[i]:=x1[i]-AA[i,j]*x0[j];
end do:
x1[i]:=x1[i]+AA[i,i]*x0[i];
x1[i]:=x1[i]/AA[i,i];
end do:
x0:=evalf(x1);
print(iter,x0);
end do:

#3 Gauss-Seidel method for matrix inverse.

> x0:=[0,0,0,0]:
for iter from 1 to 20 do
for i from 1 to n do
x1:=b[i];
for j from 1 to n do
x1:=x1-AA[i,j]*x0[j];
end do:
x1:=x1+AA[i,i]*x0[i];
x1:=x1/AA[i,i];
x0[i]:=evalf(x1);
end do:
print(iter,x0);
end do;

Geometrical meaning of Eigenvalue
固有値の
A*x = lambda*x
幾何学的意味

#1
restart;
with(LinearAlgebra):with(plots):with(plottools):

Warning, the name changecoords has been redefined

Warning, the name arrow has been redefined

> #2
A:=Matrix(1..2,1..2,[[3,2/3],[2/3,2]]);

> #3
N:=30;p1:=[];l1:=[];
for k from 0 to N-1 do
x0:=Vector([sin(2*Pi*k/N),cos(2*Pi*k/N)]);
x1:=MatrixVectorMultiply(A,x0);
p1:=[op(p1),pointplot({x0,x1})];
l1:=[op(l1),line( evalf(convert(x0,list)),evalf(convert(x1,list)) )];
end do: