国科大UCAS 2022年矩阵分析李保滨大作业matlab代码

原创，注释全，可以参考也可以拿来直接用。内含五大分解，计算矩阵行列式等matlab代码

本人使用的是MATLAB R2022b，如果你的MATLAB版本较低可能会报错。

老师给出的要求：

一个综合程序，根据选择参数的不同，实现不同的矩阵分解；

在此基础上，实现Ax=b方程组的求解，以及计算A的行列式；

Faciorization.m根据选择参数的不同，实现不同的矩阵分解。LU.m, GramSchmidt.m, Householder.m, Givens.m, URV.m分别为LU分解，Gram-Schmidt，Householder Reduction，Givens Reduction，URV相应的程序。Axb.m为求解方程Ax=b的程序。DET.m为求解方程行列式的程序。main.m为最终执行文件，具体用法均有注释。

LU.m

function[l,u,p] = LU(a)
% 示例：输入矩阵a=
%   1   2  -3   4
%   4   8  12  -8
%   2   3   2   1
%   3  -1   1  -4
[m,n]=size(a);      % 获取a的行列数

if m~=n
    error('矩阵非方阵！')
end

p = eye(n);     % 初始化P L U矩阵
l = zeros(n,n);
u = a;

for j = 1 : n-1
% 示例：第二轮循环开始时，j=2
% u =  4  8  12  −8
%      0  0  −6   6
%      0 −1  −4   5
%      0  5  10 −10
    [~, index] = max(abs(u(j:n,j)));    % 寻找当前枢轴下最大主元的相对位置索引
                                            % u(2:4,2)为[0 -1 5]'，max得到index=3
    index = index + j - 1;                  % index表示最大主元在矩阵中的行数索引
                                            % 第二轮循环中，5所在行数为4
    p([j,index],:) = p([index,j],:);        % P矩阵进行相同的行变换
                                            % p([2 4],:) = p([4 2],:)   交换p矩阵R2与R4
    u([j,index],:) = u([index,j],:);        % 行变换使当前列的最大主元到枢轴位置
                                            % u([2 4],:) = u([4 2],:)   交换u矩阵R2与R4
    l([j,index],:) = l([index,j],:);


% u =  4  8  12  −8
%      0  5  10 −10
%      0 −1  −4   5
%      0  0  −6   6

    for i = j+1 : n                         % j=2,i=3时
        multi = u(i,j)/u(j,j);              % multi = u(3,2)/u(2,2) = -1/5
        l(i,j) = multi;                     % l矩阵的(3,2)位置写入-1/5
        for k = j:n                         % 改变第i行元素，j列之前元素已经为0
            u(i,k) = u(i,k) - multi * u(j,k);
        end
% u =  4  8  12  −8
%      0  5  10 −10
%      0  0  −6   6
%      0  0  −2   3
    end

end

% 将l对角线置为1
for k = 1:n
    l(k,k)=1;
end

end

GramSchmidt.m

function[q,r] = GramSchmidt(a)

% a= 0 −20 −14
%    3  27  −4
%    4  11  −2

[m,n]=size(a);
if rank(a) ~= n
    error('!!!!!!!!!!!!矩阵列不满秩，无法分解!!!!!!!!!!!');
end

q = zeros(m,n);                     % 初始化Q R矩阵
r = zeros(n,n);

r(1,1) = norm(a(:,1));              % r(1,1) = norm([0,3,4]') = 5
q(:,1) = a(:,1)/r(1,1);             % q1=a1/5

% q=   0   0   0       r =  5  0  0
%    0.6   0   0            0  0  0
%    0.8   0   0            0  0  0

for j = 2:n                         % 从第2列到第n列
    for i = 1:j-1                   % j=2时，从第2列第1行到第2列第1行
        r(i,j) = q(:,i)' * a(:,j);  % r(i,j) = <qi|aj>
                                    % r(1,2) = -20*0+27*0.6+11*0.8 = 25
    end
    
    % 示例：j=2时结果
    % q=   0   0   0       r =  5  25 0
    %    0.6   0   0            0  0  0
    %    0.8   0   0            0  0  0

    qq = zeros(m,1);
    for i = 1: j-1
        qq = q(:,i) * r(i,j) + qq;      % qq为qi<qi|aj>对i求和
    end
    r(j,j) = norm(a(:,j)-qq);         % 对角线上元素||a2-q1<q1|a2>||
                                        % [-20,12,-9]'的模长，25
    q(:,j) = (a(:,j) - qq) / r(j,j);    % [-20,12,-9]' * 1/25 作为q的第2列

    % 示例：j=2时结果
    % q=   0   -0.8    0       r =  5  25  0
    %    0.6    0.48   0            0  25  0
    %    0.8   -0.36   0            0   0  0

end
r(abs(r)<1e-12) = 0;    % 置0

end

Householder.m

function[q,r] = Householder(a)

% 考察教材QR factorization的定义，Q为mxn维，R为nxn上三角矩阵
% 为方便计算URV分解，Householder结果为full QR factorization
% Q为mxm维，R为mxn维

% a= 0 −20 −14
%    3  27  −4
%    4  11  −2

[m,n]=size(a);
q = eye(m);
for k=1:m-1                         % 当k=2时第2轮循环
    ak = a(k:m,k:n);                % a = 5   25   −4   ak =  0  -10
                                    %     0    0  −10       -25  -10
                                    %     0  -25  −10
    mk = m-k+1;                     % mk是ak的行数
    
    e = zeros(mk, 1);
    e(1,1) = 1;                     % e=[1,0]'

    u = ak(:,1)-norm(ak(:,1))*e;    % u=25*[-1,-1]'

    R = eye(mk)-2*(u*u')/(u'*u);
    I = eye(m);
    I(k:m,k:m) = R;                 % R = 1  0
    R = I;                          %     0  R

    a = R*a;                        % …R3*R2*R1*A=T=r
    q = R*q;                        % q= R3*R2*R1

end

r=a;
q=q';   % q= (R3*R2*R1)'
r(abs(r)<1e-12) = 0;    % 置0
end

Givens.m

function[q,r] = Givens(a)

% a= 0 −20 −14
%    3  27  −4
%    4  11  −2

[m,n]=size(a);

q = eye(m); % 初始化Q矩阵

for i = 1:n
    if a(i,i)==0            % 当前位置元素为0时
        p = eye(m);
        [pivot, index] = max(abs(a(i:m,i)));    % 寻找当前枢轴下最大主元的相对位置
        index = index + i - 1;                  % index表示最大主元在矩阵中的行数索引
        p([i,index],:) = p([index,i],:);        % 将最大主元换到上面
        a = p*a;
        q = p*q;
    end
    for j = i+1:m
        c = a(i,i)/(a(i,i)^2 + a(j,i)^2)^0.5;
        s = a(j,i)/(a(i,i)^2 + a(j,i)^2)^0.5;
        p = eye(m);
                                    % p = 1  0  0
                                    %     0  1  0
                                    %     0  0  1
        p([i,j],[i,j]) = [c,s;-s,c];
                                    % p = c  s  0
                                    %    -s  c  0
                                    %     0  0  1
        a = p*a;
        q = p*q;
    end
    if i == m           % 当m<n时防止数组越界
        break;
    end
end
r = a;
q = q';
r(abs(r)<1e-12) = 0;    % 置0
end

URV.m

function[U,R,V] = URV(A)
    
    [m,n] = size(A);

    [Q1,R1] = Householder(A);   % PA = (B 0)'   其中P=Q1',[B 0]'=R1
    r = rank(R1);
    B = R1(1:r, :);

    [Q2,R2] = Householder(B');  % QB'=(T 0)' -> B=(T'|0)Q, Q=Q2'
    r = rank(R2);
    T = R2(1:r, 1:r);
    
    U = Q1;
    R = zeros(m,n);
    R(1:r,1:r) = T';
    V = Q2;

end

DET.m

function[x] = DET(A)

% det(A)=det(P)*det(L)*det(U)
% L为主对角线全为1的下三角矩阵，det(L)=1

[m,n] = size(A);
if m ~= n
    error('矩阵非方阵，无法求解！！');
end

% 进行LU分解
[~,U,P] = Factorization(A,0);


% 计算det(U)

% U为上三角矩阵,det(U)为主对角线元素之积
x = 1;
for i = 1:n
    x = x*U(i,i);
end

% 计算det(P)

% P为单位矩阵经过行变换，det(P)= 1或-1
[kk,jj] = find(P==1);   % 寻找P的元素为1的坐标，kk为所在行，jj为所在列
kk = kk';

% 计算kk的逆序数
count = 0;
for i = 1:n
    for j = i:n
        if kk(i) > kk(j)
            count = count + 1;
        end
    end
end

if mod(count,2)==1  % 逆序数为奇数时
    x = -x;
end

end

Axb.m

% 利用LU分解求Ax=b

function[x] = Axb(A,b)

[m,n] = size(A);
if m ~= n
    error('矩阵非方阵，无法求解！！');
end
if m ~= size(b,1)
    error('A与b行数不同，无法求解！！');
end

x = zeros(n,1);
y = zeros(n,1);

% 对A进行LU分解
[L,U,P] = Factorization(A,0);
b = P*b;

% 求Ly=b
% 递推公式yi=bi-(求和yj*Lij)，求和j从1到i-1
for i = 1:n
    y(i,1) = b(i,1);
    for j = 1:i-1
        y(i,1) = y(i,1)-y(j,1)*L(i,j);
    end
end

% 求解Ux=y
% 递推公式xi=(yi-(求和yj*uij))/uii，求和j从i+1到n
for i = n:-1:1  % 从n到1
    x(i,1) = y(i,1);
    for j = i+1:n
        x(i,1) = x(i,1)-x(j,1)*U(i,j);
    end
    x(i,1) = x(i,1)/U(i,i);
end

end

Factorization.m

% method为选取的分解方法
% method=0  LU分解                  [L,U,P] = Factorization(a, 0)
% method=1  GramSchmidt             [Q,R] = Factorization(a, 1)
% method=2  Householder  Reduction  [Q,R] = Factorization(a, 2)
% method=3  Givens  Reduction       [Q,R] = Factorization(a, 3)
% method=4  URV                     [U,R,V] = Factorization(a, 4)

% 考察教材QR factorization的定义，Q为mxn维，R为nxn上三角矩阵
% GramSchmidt下的QR分解采取这种形式
% 为方便进行URV分解，Householder和Givens采取full QR factorization
% 即Q为mxm维，R为mxn维


function[a1,a2,a3] = Factorization(a, method)

if method==0
    [a1,a2,a3] = LU(a);

elseif method==1
    [a1,a2] = GramSchmidt(a);

elseif method==2
    [a1,a2] = Householder(a);

elseif method==3
    [a1,a2] = Givens(a);

elseif method==4
    [a1,a2,a3] = URV(a);

else
    fprintf('输入参数有误！')

end
end

main.m

clear all
clc
%% Faciorization函数实现不同的矩阵分解
% method为选取的分解方法
% method=0  LU分解                  [L,U,P] = Factorization(a, 0)
% method=1  GramSchmidt             [Q,R] = Factorization(a, 1)
% method=2  Householder  Reduction  [Q,R] = Factorization(a, 2)
% method=3  Givens  Reduction       [Q,R] = Factorization(a, 3)
% method=4  URV                     [U,R,V] = Factorization(a, 4)

% B = [0,-20,-40;3,27,-4;4,11,-2];
B = [1,2,4,17;3,6,-12,3;2,3,-3,2;0,2,-2,6];

[L,U,P] = Factorization(B, 0);
[Q1,R1] = Factorization(B, 1);
[Q2,R2] = Factorization(B, 2);
[Q3,R3] = Factorization(B, 3);
[u,r,v] = Factorization(B, 4);

% L,U,P
% Q1,R1
% Q2,R2
% Q3,R3
% u,r,v

%% Axb函数实现求解方程Ax=b

A = [1,2,4,17;3,6,-12,3;2,3,-3,2;0,2,-2,6];
b = [17;3;3;4];

x1 = Axb(A,b);
disp(x1);

% 输出Ax=b的解x
%     2
%    -1
%     0
%     1


%% DET函数实现求矩阵行列式

A = [1,2,4,17;3,6,-12,3;2,3,-3,2;0,2,-2,6];

x2 = DET(A);
disp(x2);