一、Matlab矩阵运算
1、行列式
>> A = [1 0 2;4 7 2;6 5 7]
A =
1 0 2
4 7 2
6 5 7
>> det(A)
ans =
-5.0000
2、求逆
>> inv(A)
ans =
-7.8000 -2.0000 2.8000
3.2000 1.0000 -1.2000
4.4000 1.0000 -1.4000
3、构造增广矩阵
>> A
A =
1 0 2
4 7 2
6 5 7
>> Y = [3;4;6]
Y =
3
4
6
>> Aug = [A Y]
Aug =
1 0 2 3
4 7 2 4
6 5 7 6
4、LU分解
>> A= [1 2 6;4 8 -1;-2 3 -5]
A =
1 2 6
4 8 -1
-2 3 -5
>> [L,U,P] = lu(A)
L =
1.0000 0 0
-0.5000 1.0000 0
0.2500 0 1.0000
U =
4.0000 8.0000 -1.0000
0 7.0000 -5.5000
0 0 6.2500
P =
0 1 0
0 0 1
1 0 0
5、欧几里得范数
%行向量、列向量均可
>> Y
Y =
3
4
6
>> norm(Y)
ans =
7.8102
>> norm(Y')
ans =
7.8102
>> norm(Y,1) %1-范数
ans =
13
对norm函数的具体说明可以参照【matlab】范数小结。
二、回代法求上三角矩阵
AX=B已经进行高斯消元，得到UX=R，且对角元素非零，回代UX=R求解X：
function [ X ] = backsub( U, R )
%UX = R
%det(U)!=0
if det(U)==0
return;
end
n = length(U);
X = zeros(n,1); %列向量X
X(n) = R(n)/U(n,n);
for k=n-1:-1:1
X(k) = (R(k)-U(k,k+1:n)*X(k+1:n))/U(k,k);
end
end
测试：
>> U = [1,4,6;0,3,4;0,0,3]
U =
1 4 6
0 3 4
0 0 3
>> R = [2,4,7]
R =
2 4 7
>> backsub(U,R)
X =
0
0
2.3333
ans =
-4.8889
-1.7778
2.3333
三、高斯消元+回代
对增广矩阵[A B]变换为上三角矩阵，并回代，求出AX=B的解：
function [ X ] = uptrbk( A, B )
%AX = B
%det(A)!=0
if det(A)==0
return;
end
n = length(A);
X = zeros(n,1); %列向量X
C = zeros(1, n+1); %临时变量
Aug = [A B];
for p = 1:n-1
%偏序选主元
[Y,j] = max(abs(Aug(p:n,p)));
%行交换
C = Aug(p,:);
Aug(p,:) = Aug(j+p-1,:);
Aug(j+p-1,:) = C;
if Aug(p,p)==0
'A为奇异矩阵'
return;
end
for k = p+1:n
m = Aug(k,p)/Aug(p,p);
Aug(k,p:n+1) = Aug(k, p:n+1) - Aug(p, p:n+1)*m; %行操作
end
end
%得到上三角矩阵[U R]
X = backsub(Aug(1:n, 1:n), Aug(1:n, n+1));
%回代
end
测试：
>> A
A =
1 0 2
4 7 2
6 5 7
>> Y
Y =
3
4
6
>> uptrbk(A,Y)
ans =
-14.6000
6.4000
8.8000
四、LU分解
利用LU分解（无行交换）求解AX=B，流程如下：
- 构造LU和P，其中LU同时存储在A中，P存储在R中
- LY = PB解出Y
- UX = Y解出X
function [ X ] = lufact( A, B )
%AX = B
%det(A)!=0
if det(A)==0
return;
end
N = length(A);
X = zeros(N,1); %列向量X
Y = zeros(N,1); %列向量Y
C = zeros(1,N); %行向量C，用于行交换的中间变量
R = 1:N; %置换信息
for p = 1:N-1
%偏序选主元
[max1,j] = max(abs(A(p:N, p)));
%行交换
C = A(p,:);
A(p,:) = A(j+p-1,:);
A(j+p-1,:) = C;
%存储交换信息
d = R(p);
R(p) = R(p+j-1);
R(p+j-1) = d;
if A(p,p)==0
'A为奇异矩阵'
return;
end
for k = p+1:N;
m = A(k,p)/A(p,p);
A(k,p) = m;%利用一个矩阵A同时存储LU的信息
A(k,p+1:N) = A(k, p+1:N) - A(p, p+1:N)*m; %行操作
end
end
%得到矩阵A同时存储LU的信息
%求LY = PB
Y(1) = B(R(1));
for k = 2:N
Y(k) = B(R(k)) - A(k, 1:k-1)*Y(1:k-1);
end
%求UX = Y
X(N) = Y(N)/A(N,N);
for k=N-1:-1:1
X(k) = (Y(k)-A(k,k+1:N)*X(k+1:N))/A(k,k);
end
end
>> A
A =
1 0 2
4 7 2
6 5 7
>> Y
Y =
3
4
6
>> lufact(A, Y)
ans =
-14.6000
6.4000
8.8000
五、雅可比迭代
利用迭代求解AX = B，要求A满足严格对角优势。
function [ X ] = jacobi( A,B,P,delta,max1 )
%迭代求解AX = B，要求A满足严格对角优势
%P为X的初始值P0
%delta：允许误差范围
%max1：最大迭代次数
%det(U)!=0
N = length(A);
X = zeros(N,1);
for k = 1:max1
%计算Pk
for j = 1:N
X(j) = (B(j) - A(j, [1:j-1,j+1:N] )* P([1:j-1,j+1:N]))/A(j,j);
end
%计算误差
err = abs(norm(X - P));
P = X;
if err<delta
return;
end
end
end
>> jacobi(A,Y,[0;0;0],0.000001,10000)
ans =
0.2500
0.0833
1.4167
>> lufact(A, Y)
ans =
0.2500
0.0833
1.4167
六、高斯-赛德尔迭代
利用迭代求解AX = B，要求A满足严格对角优势。
这一迭代法收敛速度更快。
function [ X ] = gseid( A,B,P,delta,max1 )
%迭代求解AX = B，要求A满足严格对角优势
%P为X的初始值P0
%delta：允许误差范围
%max1：最大迭代次数
%det(U)!=0
N = length(A);
X = P;
for k = 1:max1
%计算Pk
for j = 1:N
X(j) = (B(j) - A(j, [1:j-1,j+1:N] )* X([1:j-1,j+1:N]))/A(j,j);
end
%计算误差
err = abs(norm(X - P));
P = X;
if err<delta
return;
end
end
>> gseid(A,Y,[0;0;0],0.000001,10000)
ans =
0.2500
0.0833
1.4167
