MATLAB中矩阵方程求解的实现

一、矩阵方程

1、定义：

 

 

2、分类

http://naotu.baidu.com/file/14d36860667d356a54490320cdab2950?token=56a299fe0e0815fe

二、M代码实现

1、M代码

function d=CDBH_for_sov_JZFC(a,b)
[m1,n1]=size(a);
[m2,n2]=size(b);
c=[a,b];
ra=rank(a);             %矩阵a的秩
rb=rank(b);             %矩阵b的秩
rc=rank(c);             %矩阵[a,b]的秩
zero=zeros(m2,n2);      %构造与b规格相同的零矩阵
pj=zero~=b;             %确定b中非零元素的个数
pj=sum(pj);
pj=sum(abs(b));         %更新，把判据改为b中所有值的绝对值之和
global i1;              %用于阶梯型的计算
global i0;              %用于阶梯型的计算，其值为当前列于当前行的差值
global cx;              %用于记录阶梯型的首个元素的位置
i1=0;
i0=0;
x_fqc=[];                  %非齐次计算中，用于记录阶梯型的首个元素的位置
x_qc=[];                  %在齐次计算中，用于记录阶梯型的首个元素的位置
cx=[];
if m1~=m2
    error('输入有误，无法计算');
    return;
end
switch pj
    %这种情况为其次方程组
    case 0
        switch rc
            %这种情况只有零解
            case n1
                
                d=zeros(n1,1);
                %disp(d);
                
                %这种情况有基础解系
            case num2cell([0:n1-1])
                %求解思路：
                %一.矩阵变换
                %1）化为行阶梯型
                %2）化为标准型
                %二.求基础解系
                %1）分别取非线性向量
                %2）则线性向量的值为，上述向量对应元素的相反数
                %三.输出结果
                
                %一、矩阵变换
                for i=1:m1-1;    %需要对行数-1行进行行变换
                    %选中了第i行
                    %如果都为非线性相关的向量，则阶梯的行列数相等，若存在线性相
                    %的向量，则需要构造一个i0来表示阶梯对应的列，并用i1表示列于
                    %行数的差值。
                    i0=i1+i;
                    %因为i1的值根据本行元素具体情况确定，因此需要注意，应当在这
                    %个for循环内完成下三角、化为1、上三角的所有操作。
                    %Step01 检查阶梯元素是否等于零
                    %       若等于零，则需要与第一个不为零的行对调，若全为零，
                    %       则i1+1，并再次循环。
                    pj_01=0;%初始化以下循环的判据
                    while pj_01==0 %利用判据pj_01来判断是否需要再次循环，在循环中，通过修改pj的值来跳出循环。
                        pj_01=1;    %没有特殊情况执行完跳出循环
                        i0=i+i1;
                        if c(i,i0)==0
                            %01 找到第i0列，i行及以下第一个非零元素的位置k
                            k=find(c([i+1:end],i0));
                            %k=k(1);
                            %02 将k行与i行对调（若k为空集，应当i1+1并再循环）
                            if k~=[];           %k不为空集时，将k行与i行对调
                                k=k(1);
                                e=c(i,:);
                                c(i,:)=c(k,:);
                                c(k,:)=e;
                                pj_01=1;
                            elseif isempty(k)==1;       %k为空集时，i1+1并再循环
                                i1=+1;
                                i0=i1+i;
                                pj_01=0;        %将判据设为0，再次循环
                            end
                        end
                    end
                    i0=i1+i;        
                    x_qc=[x_qc,i0];     %将此时的列位置进行记录
                    %Step2 检查阶梯元素是否为1，若不为1，则将其化为1
                    if c(i,i0)~=0&&c(i,i0)~=0;
                        c(i,:)=c(i,:)/c(i,i0);
                    end
                    %Step3 将i0列，第i行一下的元素消减为0
                    for j=i+1:m1
                        if c(j,i0)~=0
                            c(j,:)=c(j,:)-c(j,i0)*c(i,:);
                        end
                    end
                    %Step4 将i0列，第i行一上的元素消减为0
                    for j=1:i-1
                        if c(j,i0)~=0
                            c(j,:)=c(j,:)-c(j,i0)*c(i,:);
                        end
                    end
                end
                %更新！检查最后一行
                if sum(abs(c(m1,[1:n1])))~=0
                    c(m1,:)=c(m1,:)/c(m1,n1);
                    for j=1:m1-1
                        if c(j,n1)~=0
                            c(j,:)=c(j,:)-c(j,n1)*c(m1,:);
                        end
                    end
                end
                
                %成功化为标准型
                %disp(c)    
                %二、求基础解系
                %依次选定线性相关向量所在列，查找对应非线性相关的元素的值，并将其取负，放入解系矩阵。
                %Step1 根据x_qc构造线性相关向量位置向量，构造基础解系矩阵
                no_x_qc=ones(1,n1);     %初始化与a列咧数相等的1向量
                no_x_qc(x_qc)=0;        %其中线性无关向量位置设为0
                no_x_qc=find(no_x_qc);  %找到非零元素位置，命名为no_x_qc
                jcjx=zeros(n1,rc);      %初始化基础解系（行：A的列，列：C的秩）
                %Step2 选定线性相关所在列（使用for循环）
                for i=1:length(no_x_qc)
                    
                    %Step3 查找该列中，线性无关向量对应的元素的值
                    for j=1:length(x_qc)
                        Psi=c(j,no_x_qc(i));
                        %Step4 将查找到的值取负，并放入基础解系的第i列的相应位置
                        Psi=Psi*-1;
                        jcjx(x_qc(j),i)=Psi;
                    end
                    %Step5 将基础解系中，当前列对应的位置的元素赋值为1
                    jcjx(no_x_qc(i),i)=1;
                end
                %disp(jcjx)
                %三、输出结果
                disp '齐次型，结果为基础解系'
                d=jcjx;
                disp 'x1'
                disp '... = k1*psi1+...+kr*psir'
                disp 'xn'
                    otherwise
                        error('输入有误，无法计算');
                        return;
                end
                %这种话情况为非齐次方程
    otherwise
        %这种情况无解
        if ra<rc
            error('非齐次方程组无解');
            return;
            %这种情况唯一解
        elseif ra==rc&&rc==n1   %这时，a必然为方阵，但是可以使用上述齐次方程的解法
            %思路：
            %一、矩阵变换
            %二、求解
            %三、输出结果
            
            %一、矩阵变换
            for i=1:m1-1;    %需要对行数-1行进行行变换
                %选中了第i行
                %如果都为非线性相关的向量，则阶梯的行列数相等，若存在线性相
                %的向量，则需要构造一个i0来表示阶梯对应的列，并用i1表示列于
                %行数的差值。
                i0=i1+i;
                %因为i1的值根据本行元素具体情况确定，因此需要注意，应当在这
                %个for循环内完成下三角、化为1、上三角的所有操作。
                %Step01 检查阶梯元素是否等于零
                %       若等于零，则需要与第一个不为零的行对调，若全为零，
                %       则i1+1，并再次循环。
                pj_01=0;%初始化以下循环的判据
                while pj_01==0 %利用判据pj_01来判断是否需要再次循环，在循环中，通过修改pj的值来跳出循环。
                    pj_01=1;    %没有特殊情况执行完跳出循环
                    i0=i+i1;
                    if c(i,i0)==0
                        %01 找到第i0列，i行及以下第一个非零元素的位置k
                        k=find(c([i+1:end],i0));
                        %k=k(1);
                        %02 将k行与i行对调（若k为空集，应当i1+1并再循环）
                        if k~=[];           %k不为空集时，将k行与i行对调
                            k=k(1);
                            e=c(i,:);
                            c(i,:)=c(k,:);
                            c(k,:)=e;
                            pj_01=1;
                        elseif isempty(k)==1;       %k为空集时，i1+1并再循环
                            i1=+1;
                            i0=i1+i;
                            pj_01=0;        %将判据设为0，再次循环
                        end
                    end
                end
                i0=i1+i;
                x_qc=[x_qc,i0];     %将此时的列位置进行记录
                %Step2 检查阶梯元素是否为1，若不为1，则将其化为1
                if c(i,i0)~=0&&c(i,i0)~=0;
                    c(i,:)=c(i,:)/c(i,i0);
                end
                %Step3 将i0列，第i行一下的元素消减为0
                for j=i+1:m1
                    if c(j,i0)~=0
                        c(j,:)=c(j,:)-c(j,i0)*c(i,:);
                    end
                end
                %Step4 将i0列，第i行一上的元素消减为0
                for j=1:i-1
                    if c(j,i0)~=0
                        c(j,:)=c(j,:)-c(j,i0)*c(i,:);
                    end
                end
            end
            %更新！检查最后一行
            if sum(abs(c(m1,[1:n1])))~=0
                c(m1,:)=c(m1,:)/c(m1,n1);
                for j=1:m1-1
                    if c(j,n1)~=0
                        c(j,:)=c(j,:)-c(j,n1)*c(m1,:);
                    end
                end
            end
            
            %成功化为标准型
            %disp(c)
            %二、求解
            %c变换后，其b的位置的数据即为解
            jx=c(:,[n1+1:end]);
            %三、输出结果
            d=jx;
            disp '非齐次型，唯一解'
            
            %无穷多解
        elseif ra==rc&&rc<n1
            %思路：
            %一、矩阵变换（利用齐次方程的代码）
            for i=1:m1-1;    %需要对行数-1行进行行变换
                %选中了第i行
                %如果都为非线性相关的向量，则阶梯的行列数相等，若存在线性相
                %的向量，则需要构造一个i0来表示阶梯对应的列，并用i1表示列于
                %行数的差值。
                i0=i1+i;
                %因为i1的值根据本行元素具体情况确定，因此需要注意，应当在这
                %个for循环内完成下三角、化为1、上三角的所有操作。
                %Step01 检查阶梯元素是否等于零
                %       若等于零，则需要与第一个不为零的行对调，若全为零，
                %       则i1+1，并再次循环。
                pj_01=0;%初始化以下循环的判据
                while pj_01==0 %利用判据pj_01来判断是否需要再次循环，在循环中，通过修改pj的值来跳出循环。
                    pj_01=1;    %没有特殊情况执行完跳出循环
                    i0=i+i1;
                    if c(i,i0)==0
                        %01 找到第i0列，i行及以下第一个非零元素的位置k
                        k=find(c([i+1:end],i0));
                        %k=k(1);
                        %02 将k行与i行对调（若k为空集，应当i1+1并再循环）
                        if k~=[];           %k不为空集时，将k行与i行对调
                            k=k(1);
                            e=c(i,:);
                            c(i,:)=c(k,:);
                            c(k,:)=e;
                            pj_01=1;
                        elseif isempty(k)==1;       %k为空集时，i1+1并再循环
                            i1=+1;
                            i0=i1+i;
                            pj_01=0;        %将判据设为0，再次循环
                        end
                    end
                end
                i0=i1+i;
                x_qc=[x_qc,i0];     %将此时的列位置进行记录
                %Step2 检查阶梯元素是否为1，若不为1，则将其化为1
                if c(i,i0)~=0&&c(i,i0)~=0;
                    c(i,:)=c(i,:)/c(i,i0);
                end
                %Step3 将i0列，第i行一下的元素消减为0
                for j=i+1:m1
                    if c(j,i0)~=0
                        c(j,:)=c(j,:)-c(j,i0)*c(i,:);
                    end
                end
                %Step4 将i0列，第i行一上的元素消减为0
                for j=1:i-1
                    if c(j,i0)~=0
                        c(j,:)=c(j,:)-c(j,i0)*c(i,:);
                    end
                end
            end
            %更新！检查最后一行
            if sum(abs(c(m1,[1:n1])))~=0
                c(m1,:)=c(m1,:)/c(m1,n1);
                for j=1:m1-1
                    if c(j,n1)~=0
                        c(j,:)=c(j,:)-c(j,n1)*c(m1,:);
                    end
                end
            end
            
            %成功化为标准型
            %disp(c)
            %二、求特解Psit
            x_fqc=x_qc;                 %为了便于阅读，使用x_fqc代替x_qc
            Psit=zeros(1,n1);           %初始化特解（长度：a的列数）
            nx=length(x_fqc);           %nx为线性无关向量的数量
            for i=1:nx
                ip=x_fqc(i);
                Psit1=c(:,[n1+1:end]);
                Psit(ip)=Psit1(i);
                disp(Psit);
            end
            Psit=Psit';
            %三、构造对应齐次方程，并解出基础解系
            b1=zeros(size(b));
            
            px=ones(1,n1);
            px(x_fqc)=0;
            px=find(px);
            npx=length(x_fqc);
            a1=a;
            %for i=1:npx
            %    ipx=px(i);
            %    a1(:,ipx)=a1(:,ipx)*-1;
            %    disp(a1);
            %end
            
            d1=CDBH_for_sov_JZFC(a1,b1);
            %四、构造通解
            tj=[d1,Psit];
            %五、输出结果
            d=tj;
            disp '非齐次，一般解：'
            disp 'x=k1*psi1+...+kn*psin+psit'
        else
            error('输入有误，无法计算');
            return;
        end
end       

 

2、例子（用于验证）

（1）齐次方程（无限解）：

a=[1,1,-1,-1;2,-5,3,2;7,-7,3,1];
b=[0,0,0]';

结果：

 

（2）非齐次方程（唯一解）：

a=[89,14,81,19;95,25,24,25;54,84,92,61;13,25,34,47];
b=[436,317,742,353]';

结果：

>>x =1       2       3       4 

 

（3）非齐次方程（无限解）：

a=[1,-1,-1,1;1,-1,1,-3;1,-1,-2,3];
b=[0,1,-1/2]';

结果：