Maple,编写一个 Proc 来计算矩阵的逆
Maple, writing a Proc to caculate inverse of a Matrix
所以我正在研究一个由两部分组成的问题。我确实在 :this useful Forum 的帮助下完成了第一部分。有些人已经尝试解决问题的第一部分,我接受了他们的代码。
问题:
- 在 Maple 中编写一个过程(我在这里命名为 "reduced"),计算矩阵的简化阶梯形式。
- 编写一个使用“约简”计算矩阵逆的过程。
第一部分的代码(代码已经过测试,我声称它可以正确运行)
multiplikation:= proc(
m::posint,
a::depends(And(posint, satisfies(a-> a <= m))),
b::And({float, rational}, Not(identical(0,0.)))
)
Matrix((m,m), (i,j)-> `if`(i=j, `if`(i=a, b, 1), 0))
end proc:
addition:= proc(
m::posint,
a::depends(And(posint, satisfies(a-> a <= m))),
b::depends(And(posint, satisfies(b-> b <= m))),
c::And({float, rational}, Not(identical(0,0.)))
)
Matrix((m,m), (i,j)-> `if`(i=a and j=b, c, `if`(i=j, 1, 0)))
end proc:
perm:= proc(
m::posint,
a::depends(And(posint, satisfies(a-> a <= m))),
b::depends(And(posint, satisfies(b-> b <= m and a<>b)))
)
Matrix((m,m), (i,j)-> `if`({i,j}={a,b} or i=j and not i in {a,b}, 1, 0))
end proc:
和主要过程:
reduced:= proc(B::Matrix)
uses LA= LinearAlgebra;
local
M:= B, l:= 1, #l is current column.
m:= LA:-RowDimension(M), n:= LA:-ColumnDimension(M), i, j
;
for i to m do #going through every row item
#l needs to be less than column number n.
if n < l then return M end if;
j:= i; #Initialize current row number.
while M[j,l]=0 do #Search for 1st row item <> 0.
j:= j+1;
if m < j then #End of row: Go to next column.
j:= i;
l:= l+1;
if n < l then return M fi #end of column and row
end if
end do;
if j<>i then M:= perm(m,j,i).M end if; #Permute rows j and i
#Multiply row i with 1/M[i,l], if it's not 0.
if M[i,l] <> 0 then M:= multiplikation(m,i,1/M[i,l]).M fi;
#Subtract each row j with row i for M[j,l]-times.
for j to m do if j<>i then M:= addition(m,j,i,-M[j,l]).M fi od;
l:= l+1 #Increase l by 1; next iteration i increase either.
end do;
return M
end proc:
如果您需要有关上述代码的任何其他信息,我会解释更多。
对于第二部分,我想使用 Gauss Jordan Algorithmus,但我有一个问题:
我不能将单位矩阵用作“减少”中的参数。因为它的行数和列数都是 0。
你知道我如何在我的 proc:reduced 的帮助下实现 Gauss Jordan 算法吗?
既定目标是利用 reduced
程序。
一种方法是通过单位矩阵扩充输入矩阵,减少它,然后 return 扩充矩阵的右半部分。
将输入矩阵转换为单位矩阵的步骤ps也将单位矩阵转换为(输入矩阵的)逆矩阵。
例如,使用您的程序,
inv := proc(B::Matrix(square))
local augmented,m;
uses LinearAlgebra;
m := RowDimension(B);
augmented := <<B|IdentityMatrix(m)>>;
return reduced(augmented)[..,m+1..-1];
end proc:
MM := LinearAlgebra:-RandomMatrix(3,generator=1..5);
[1 4 2]
[ ]
MM := [1 5 3]
[ ]
[2 3 5]
ans := inv(MM);
[ 8 -7 1]
[ - -- -]
[ 3 3 3]
[ ]
[ 1 1 -1]
ans := [ - - --]
[ 6 6 6 ]
[ ]
[-7 5 1]
[-- - -]
[6 6 6]
ans.MM, MM.ans;
[1 0 0] [1 0 0]
[ ] [ ]
[0 1 0], [0 1 0]
[ ] [ ]
[0 0 1] [0 0 1]
ps。您可能还想考虑当矩阵不可逆时您的 reduce
程序会做什么。
所以我正在研究一个由两部分组成的问题。我确实在 :this useful Forum 的帮助下完成了第一部分。有些人已经尝试解决问题的第一部分,我接受了他们的代码。
问题:
- 在 Maple 中编写一个过程(我在这里命名为 "reduced"),计算矩阵的简化阶梯形式。
- 编写一个使用“约简”计算矩阵逆的过程。
第一部分的代码(代码已经过测试,我声称它可以正确运行)
multiplikation:= proc(
m::posint,
a::depends(And(posint, satisfies(a-> a <= m))),
b::And({float, rational}, Not(identical(0,0.)))
)
Matrix((m,m), (i,j)-> `if`(i=j, `if`(i=a, b, 1), 0))
end proc:
addition:= proc(
m::posint,
a::depends(And(posint, satisfies(a-> a <= m))),
b::depends(And(posint, satisfies(b-> b <= m))),
c::And({float, rational}, Not(identical(0,0.)))
)
Matrix((m,m), (i,j)-> `if`(i=a and j=b, c, `if`(i=j, 1, 0)))
end proc:
perm:= proc(
m::posint,
a::depends(And(posint, satisfies(a-> a <= m))),
b::depends(And(posint, satisfies(b-> b <= m and a<>b)))
)
Matrix((m,m), (i,j)-> `if`({i,j}={a,b} or i=j and not i in {a,b}, 1, 0))
end proc:
和主要过程:
reduced:= proc(B::Matrix)
uses LA= LinearAlgebra;
local
M:= B, l:= 1, #l is current column.
m:= LA:-RowDimension(M), n:= LA:-ColumnDimension(M), i, j
;
for i to m do #going through every row item
#l needs to be less than column number n.
if n < l then return M end if;
j:= i; #Initialize current row number.
while M[j,l]=0 do #Search for 1st row item <> 0.
j:= j+1;
if m < j then #End of row: Go to next column.
j:= i;
l:= l+1;
if n < l then return M fi #end of column and row
end if
end do;
if j<>i then M:= perm(m,j,i).M end if; #Permute rows j and i
#Multiply row i with 1/M[i,l], if it's not 0.
if M[i,l] <> 0 then M:= multiplikation(m,i,1/M[i,l]).M fi;
#Subtract each row j with row i for M[j,l]-times.
for j to m do if j<>i then M:= addition(m,j,i,-M[j,l]).M fi od;
l:= l+1 #Increase l by 1; next iteration i increase either.
end do;
return M
end proc:
如果您需要有关上述代码的任何其他信息,我会解释更多。
对于第二部分,我想使用 Gauss Jordan Algorithmus,但我有一个问题: 我不能将单位矩阵用作“减少”中的参数。因为它的行数和列数都是 0。
你知道我如何在我的 proc:reduced 的帮助下实现 Gauss Jordan 算法吗?
既定目标是利用 reduced
程序。
一种方法是通过单位矩阵扩充输入矩阵,减少它,然后 return 扩充矩阵的右半部分。
将输入矩阵转换为单位矩阵的步骤ps也将单位矩阵转换为(输入矩阵的)逆矩阵。
例如,使用您的程序,
inv := proc(B::Matrix(square))
local augmented,m;
uses LinearAlgebra;
m := RowDimension(B);
augmented := <<B|IdentityMatrix(m)>>;
return reduced(augmented)[..,m+1..-1];
end proc:
MM := LinearAlgebra:-RandomMatrix(3,generator=1..5);
[1 4 2]
[ ]
MM := [1 5 3]
[ ]
[2 3 5]
ans := inv(MM);
[ 8 -7 1]
[ - -- -]
[ 3 3 3]
[ ]
[ 1 1 -1]
ans := [ - - --]
[ 6 6 6 ]
[ ]
[-7 5 1]
[-- - -]
[6 6 6]
ans.MM, MM.ans;
[1 0 0] [1 0 0]
[ ] [ ]
[0 1 0], [0 1 0]
[ ] [ ]
[0 0 1] [0 0 1]
ps。您可能还想考虑当矩阵不可逆时您的 reduce
程序会做什么。