MATLAB基础篇——线性代数应用

在MATLAB基础篇——基础语法介绍了有关矩阵的基本运算和线性代数中的一些基本问题。这里我们再讨论如下问题：1.求向量组的极大线性无关组；2.化二次型为标准形；3.判断二次型的正定性；4.解线性方程组

向量组的极大线性无关组

通过向量组的秩来讨论向量组的线性相关性。由于矩阵的秩=行秩=列秩，所以求向量组的秩相当于求由向量组作为列（行）向量构建的矩阵的秩——通过rank函数，再利用rref函数可以得到该矩阵的标准阶梯形矩阵，从而得到向量组的极大线性无关组，及其他向量用该极大线性无关组的表示

a1=[2 1 3 2];
a2=[3 2 -2 -3];
a3=[1 0 8 7];
a4=[-3 -2 3 4];
a5=[-7 -4 0 3];
A=[a1;a2;a3;a4;a5]';  %转置作为列向量
rank(A)
rref(A)   

% 秩为3
ans =

     3

% a1' ,a2',a4'作为极大线性无关组，a3=2*a1-a2,a5=-2*a1+3*a2+4*a4
ans =

     1     0     2     0    -2
     0     1    -1     0     3
     0     0     0     1     4
     0     0     0     0     0

化二次型为标准形

线代里面，我们学习到化二次型$f=x^{T}Ax$为标准形$f=y^{T}By$的过程就是二次型的实对称矩阵$A$正交相似化的过程，二次型标准形的系数就是$A$的特征值${\lambda }_1,...,{\lambda }_n$，即$B=diag({\lambda }_1,...,{\lambda }_n)$,由$A$的特征值构成的对角阵，而该线性变换$x=Py$中的$P$为特征值对应的特征向量，所以该问题转化为求实对称矩阵$A$的特征值和特征向量——通过eig函数

%eig(A),e=eig(A) 得到A的特征值构成的列向量
%[P,E]=eig(A) P为特征向量作为列向量构成的矩阵，E为对角线元素为特征值的对角阵
%注意，若A为实对称矩阵，则得到的特征向量矩阵是正交矩阵

clear
clc
A=[-5 2 2;2 -6 0;2 0 -4];  %二次型的实对称矩阵
[P,E]=eig(A);
syms y1 y2 y3
y=[y1;y2;y3];
x=P*y; %线性变换（正交变换）
f1=x'*A*x; %原二次型
f2=y'*E*y;  %标准二次型

二次型的正定性

在线性代数的学习中，可以知道二次型的正定性判断可以通过二次型的实对称矩阵的特征值或者顺序主子式来判断。

二次型的分类：

正定性	定义	特征值	顺序主子式
正定	$\forall$ x≠0，f=$x^{T}Ax>0$	$\forall {\lambda }_k$＞0	$\forall A_k$>0
负定	$\forall$ x≠0，f=$x^{T}Ax<0$	$\forall {\lambda }_k$<0	$\forall A_k$>0,$k$为偶， $\forall A_k$<0,$k$为奇
半正定	$\forall$ x≠0，f=$x^{T}Ax\ge 0$	$\forall {\lambda }_k$≥0	$\forall A_k$≥0
半负定	$\forall$ x≠0，f=$x^{T}Ax\le 0$	$\forall {\lambda }_k$≤0	$\forall A_k$≥0,$k$为偶， $\forall A_k$≤0,$k$为奇
不定	$\forall$ x≠0,存在f>0,也存在f<0	<0，>0	不满足其他类型即是不定

clear
clc
A=[2 0 0;0 3 2;0 2 3];
%特征值法
e=eig(A)
if all(e>0)
    fprintf('正定');
elseif all(e<0)
    fprintf('负定');
elseif all(e>=0)
    fprintf('半正定');
elseif all(e<=0)
    fprintf('半负定');
else
    fprintf('不定');
end

%顺序主子式
Ak=[];
for i=1:length(A)
    Ak(i)=det(A(1:i,1:i));
end
Ak  %   2     6    10都大于0，为正定


%改进版
clear
n=length(A);
d=zeros(1,n);
for i=1:n
    d(i)=det(A(1:i,1:i));
end
d
c1=zeros(1,3);%记录偶顺序主子式>0,<0,=0
c2=zeros(1,3);%记录奇顺序主子式
for i=1:n
    if mod(i,2)==0
        if d(i)>0
            c1(1)=c1(1)+1;
        elseif d(i)<0
            c1(2)=c1(2)+1;
        else
            c1(3)=c1(3)+1;
        end
    else
         if d(i)>0
            c2(1)=c2(1)+1;
        elseif d(i)<0
            c2(2)=c2(2)+1;
        else
            c2(3)=c2(3)+1;
        end
    end   
end
c1
c2
if c1(1)+c2(1)==n
    disp('正定')
elseif c1(1)==fix(n/2) & c2(2)==n-fix(n/2)
    disp('负定')
elseif c1(1)+c1(3)+c2(1)+c2(3)==n
    disp('半正定')
elseif c1(1)+c1(3)==fix(n/2) & c2(2)+c2(3)==n-fix(n/2)
    disp('半负定')
else
    disp('不定')
end

解线性方程组

在线性代数的学习中，我们知道方程组解的情况及其判断，设$A$为线性方程组的系数矩阵，$B$为线性方程组的增广矩阵，方程组变量个数为n（讨论的是方程个数$\le n$）则对于方程组$Ax=b$
$case1：$当$R(A)=R(B)=n$，有唯一解,
$case2：$当$R(A)=R(B)<n$，有无穷多个解,
$case3：$当$R(A)\ne R(B)$，无解
则求线性方程组解的问题转化为求方程组系数矩阵及增广矩阵的秩的问题
对于$case1,A$满秩，直接利用$x=A$ \ $b$或$inv(A)\ast b$ 得到解
对于$case2$, 可以利用$rref(B)$得到阶梯形方程组的增广矩阵，进一步得到方程的基础解系

clear
clc
A=[1 -2 3 1 1;1 1 -1 -1 -2;2 -1 1 0 -2;2 2 5 -1 1];
b=[7;2;7;18];
B=[A b];
n=length(A);
RA=rank(A);RB=rank(B);
if RA==RB
    if RA==n
        fprintf('有唯一解')
        x=A\b
    else
        fprintf('有无穷多个解')
        d=rref(B)
    end
else
    fprintf('无解');
end


% 求解结果
有无穷多个解
d =

     1     0     0     0    -2     3
     0     1     0     0    -1     1
     0     0     1     0     1     2
     0     0     0     1    -2     0

当一个齐次线性方程组有无穷多个解时，可以利用null函数求其基础解系，null返回的矩阵的列向量组即为方程组的基础解系，而有无穷多个解的非齐次方程组可以利用pinv函数求得一个特解，再利用null求相应齐次方程组的基础解系，从而得到非齐次方程组的通解

%齐次方程组
clear
clc
A=[1 1 1 4 -3;2 1 3 5 -5;1 -1 3 -2 -1;3 1 5 6 -7]; %齐次方程组系数矩阵
if rank(A)==length(A)
    fprintf('有唯一解，即零解')
else
    fprintf('有无穷多个解，即有非零解')
    K=null(A,'r')
    %把解表示成x=k_1x_1+k_2x_2+...+k_nx_n
    syms k 
    x=0;
    for i=1:length(A)-rank(A)
        k(i)=strcat('k',num2str(i));
        x=x+k(i)*K(:,i);
    end
    x
end

%结果
x =
 
 2*k3 - k2 - 2*k1
   k1 - 3*k2 + k3
               k1
               k2
               k3

% 非齐次方程组
clear
clc
A=[1 5 -1 -1;1 -2 1 3;3 8 -1 1;1 -9 3 7];
b=[-1;3;1;7];
B=[A b];
RA=rank(A);
RB=rank(B);
n=length(A(1,:));
if RA==RB
    if RA==n
        fprintf('有唯一解')
        x=A\b;
    else
        fprintf('有无穷多个解')
        r=pinv(A)*b;  %特解
        K=null(A,'r'); %齐次方程组的基础解系
        %把解表示成x=r+k_1x_1+k_2x_2+...+k_nx_n
        syms k 
        x=r;
        for i=1:length(K(1,:))
            k(i)=strcat('k',num2str(i));
            x=x+k(i)*K(:,i);
        end
        x
    end
else
    fprintf('无解')
end