MATLAB中矩阵方程求解的实现
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]';