matlab特征值是空集,MATLAB中矩阵方程求解的实现

1 function d=CDBH_for_sov_JZFC(a,b)2 [m1,n1]=size(a);3 [m2,n2]=size(b);4 c=[a,b];5 ra=rank(a); %矩阵a的秩6 rb=rank(b); %矩阵b的秩7 rc=rank(c); %矩阵[a,b]的秩8 zero=zeros(m2,n2); %构造与b规格相同的零矩阵9 pj=zero~=b; %确定b中非零元素的个数10 pj=sum(pj);11 pj=sum(abs(b)); %更新，把判据改为b中所有值的绝对值之和12 global i1; %用于阶梯型的计算13 global i0; %用于阶梯型的计算，其值为当前列于当前行的差值14 global cx; %用于记录阶梯型的首个元素的位置15 i1=0;16 i0=0;17 x_fqc=[]; %非齐次计算中，用于记录阶梯型的首个元素的位置18 x_qc=[]; %在齐次计算中，用于记录阶梯型的首个元素的位置19 cx=[];20 if m1~=m221 error('输入有误，无法计算');22 return;23 end24 switchpj25 %这种情况为其次方程组26 case 0

27 switchrc28 %这种情况只有零解29 casen130

31 d=zeros(n1,1);32 %disp(d);33

34 %这种情况有基础解系35 case num2cell([0:n1-1])36 %求解思路：37 %一.矩阵变换38 %1)化为行阶梯型39 %2)化为标准型40 %二.求基础解系41 %1)分别取非线性向量42 %2)则线性向量的值为，上述向量对应元素的相反数43 %三.输出结果44

45 %一、矩阵变换46 for i=1:m1-1; %需要对行数-1行进行行变换47 %选中了第i行48 %如果都为非线性相关的向量，则阶梯的行列数相等，若存在线性相49 %的向量，则需要构造一个i0来表示阶梯对应的列，并用i1表示列于50 %行数的差值。51 i0=i1+i;52 %因为i1的值根据本行元素具体情况确定，因此需要注意，应当在这53 %个for循环内完成下三角、化为1、上三角的所有操作。54 %Step01 检查阶梯元素是否等于零55 %若等于零，则需要与第一个不为零的行对调，若全为零，56 % 则i1+1，并再次循环。57 pj_01=0;%初始化以下循环的判据58 while pj_01==0 %利用判据pj_01来判断是否需要再次循环，在循环中，通过修改pj的值来跳出循环。59 pj_01=1; %没有特殊情况执行完跳出循环60 i0=i+i1;61 if c(i,i0)==0

62 %01找到第i0列，i行及以下第一个非零元素的位置k63 k=find(c([i+1:end],i0));64 %k=k(1);65 %02 将k行与i行对调(若k为空集，应当i1+1并再循环)66 if k~=[]; %k不为空集时，将k行与i行对调67 k=k(1);68 e=c(i,:);69 c(i,:)=c(k,:);70 c(k,:)=e;71 pj_01=1;72 elseif isempty(k)==1; %k为空集时，i1+1并再循环73 i1=+1;74 i0=i1+i;75 pj_01=0; %将判据设为0，再次循环76 end77 end78 end79 i0=i1+i;80 x_qc=[x_qc,i0]; %将此时的列位置进行记录81 %Step2 检查阶梯元素是否为1，若不为1，则将其化为182 if c(i,i0)~=0&&c(i,i0)~=0;83 c(i,:)=c(i,:)/c(i,i0);84 end85 %Step3 将i0列，第i行一下的元素消减为086 for j=i+1:m187 if c(j,i0)~=0

88 c(j,:)=c(j,:)-c(j,i0)*c(i,:);89 end90 end91 %Step4 将i0列，第i行一上的元素消减为092 for j=1:i-1

93 if c(j,i0)~=0

94 c(j,:)=c(j,:)-c(j,i0)*c(i,:);95 end96 end97 end98 %更新！检查最后一行99 if sum(abs(c(m1,[1:n1])))~=0

100 c(m1,:)=c(m1,:)/c(m1,n1);101 for j=1:m1-1

102 if c(j,n1)~=0

103 c(j,:)=c(j,:)-c(j,n1)*c(m1,:);104 end105 end106 end107

108 %成功化为标准型109 %disp(c)110 %二、求基础解系111 %依次选定线性相关向量所在列，查找对应非线性相关的元素的值，并将其取负，放入解系矩阵。112 %Step1 根据x_qc构造线性相关向量位置向量，构造基础解系矩阵113 no_x_qc=ones(1,n1); %初始化与a列咧数相等的1向量114 no_x_qc(x_qc)=0; %其中线性无关向量位置设为0115 no_x_qc=find(no_x_qc); %找到非零元素位置，命名为no_x_qc116 jcjx=zeros(n1,rc); %初始化基础解系(行：A的列，列：C的秩)117 %Step2 选定线性相关所在列(使用for循环)118 for i=1:length(no_x_qc)119

120 %Step3 查找该列中，线性无关向量对应的元素的值121 for j=1:length(x_qc)122 Psi=c(j,no_x_qc(i));123 %Step4 将查找到的值取负，并放入基础解系的第i列的相应位置124 Psi=Psi*-1;125 jcjx(x_qc(j),i)=Psi;126 end127 %Step5 将基础解系中，当前列对应的位置的元素赋值为1128 jcjx(no_x_qc(i),i)=1;129 end130 %disp(jcjx)131 %三、输出结果132 disp '齐次型，结果为基础解系'

133 d=jcjx;134 disp 'x1'

135 disp '... = k1*psi1+...+kr*psir'

136 disp 'xn'

137 otherwise138 error('输入有误，无法计算');139 return;140 end141 %这种话情况为非齐次方程142 otherwise143 %这种情况无解144 if ra

154 %一、矩阵变换155 for i=1:m1-1; %需要对行数-1行进行行变换156 %选中了第i行157 %如果都为非线性相关的向量，则阶梯的行列数相等，若存在线性相158 %的向量，则需要构造一个i0来表示阶梯对应的列，并用i1表示列于159 %行数的差值。160 i0=i1+i;161 %因为i1的值根据本行元素具体情况确定，因此需要注意，应当在这162 %个for循环内完成下三角、化为1、上三角的所有操作。163 %Step01 检查阶梯元素是否等于零164 %若等于零，则需要与第一个不为零的行对调，若全为零，165 % 则i1+1，并再次循环。166 pj_01=0;%初始化以下循环的判据167 while pj_01==0 %利用判据pj_01来判断是否需要再次循环，在循环中，通过修改pj的值来跳出循环。168 pj_01=1; %没有特殊情况执行完跳出循环169 i0=i+i1;170 if c(i,i0)==0

171 %01找到第i0列，i行及以下第一个非零元素的位置k172 k=find(c([i+1:end],i0));173 %k=k(1);174 %02 将k行与i行对调(若k为空集，应当i1+1并再循环)175 if k~=[]; %k不为空集时，将k行与i行对调176 k=k(1);177 e=c(i,:);178 c(i,:)=c(k,:);179 c(k,:)=e;180 pj_01=1;181 elseif isempty(k)==1; %k为空集时，i1+1并再循环182 i1=+1;183 i0=i1+i;184 pj_01=0; %将判据设为0，再次循环185 end186 end187 end188 i0=i1+i;189 x_qc=[x_qc,i0]; %将此时的列位置进行记录190 %Step2 检查阶梯元素是否为1，若不为1，则将其化为1191 if c(i,i0)~=0&&c(i,i0)~=0;192 c(i,:)=c(i,:)/c(i,i0);193 end194 %Step3 将i0列，第i行一下的元素消减为0195 for j=i+1:m1196 if c(j,i0)~=0

197 c(j,:)=c(j,:)-c(j,i0)*c(i,:);198 end199 end200 %Step4 将i0列，第i行一上的元素消减为0201 for j=1:i-1

202 if c(j,i0)~=0

203 c(j,:)=c(j,:)-c(j,i0)*c(i,:);204 end205 end206 end207 %更新！检查最后一行208 if sum(abs(c(m1,[1:n1])))~=0

209 c(m1,:)=c(m1,:)/c(m1,n1);210 for j=1:m1-1

211 if c(j,n1)~=0

212 c(j,:)=c(j,:)-c(j,n1)*c(m1,:);213 end214 end215 end216

217 %成功化为标准型218 %disp(c)219 %二、求解220 %c变换后，其b的位置的数据即为解221 jx=c(:,[n1+1:end]);222 %三、输出结果223 d=jx;224 disp '非齐次型，唯一解'

225

226 %无穷多解227 elseif ra==rc&&rc

246 %01找到第i0列，i行及以下第一个非零元素的位置k247 k=find(c([i+1:end],i0));248 %k=k(1);249 %02 将k行与i行对调(若k为空集，应当i1+1并再循环)250 if k~=[]; %k不为空集时，将k行与i行对调251 k=k(1);252 e=c(i,:);253 c(i,:)=c(k,:);254 c(k,:)=e;255 pj_01=1;256 elseif isempty(k)==1; %k为空集时，i1+1并再循环257 i1=+1;258 i0=i1+i;259 pj_01=0; %将判据设为0，再次循环260 end261 end262 end263 i0=i1+i;264 x_qc=[x_qc,i0]; %将此时的列位置进行记录265 %Step2 检查阶梯元素是否为1，若不为1，则将其化为1266 if c(i,i0)~=0&&c(i,i0)~=0;267 c(i,:)=c(i,:)/c(i,i0);268 end269 %Step3 将i0列，第i行一下的元素消减为0270 for j=i+1:m1271 if c(j,i0)~=0

272 c(j,:)=c(j,:)-c(j,i0)*c(i,:);273 end274 end275 %Step4 将i0列，第i行一上的元素消减为0276 for j=1:i-1

277 if c(j,i0)~=0

278 c(j,:)=c(j,:)-c(j,i0)*c(i,:);279 end280 end281 end282 %更新！检查最后一行283 if sum(abs(c(m1,[1:n1])))~=0

284 c(m1,:)=c(m1,:)/c(m1,n1);285 for j=1:m1-1

286 if c(j,n1)~=0

287 c(j,:)=c(j,:)-c(j,n1)*c(m1,:);288 end289 end290 end291

292 %成功化为标准型293 %disp(c)294 %二、求特解Psit295 x_fqc=x_qc; %为了便于阅读，使用x_fqc代替x_qc296 Psit=zeros(1,n1); %初始化特解(长度：a的列数)297 nx=length(x_fqc); %nx为线性无关向量的数量298 for i=1:nx299 ip=x_fqc(i);300 Psit1=c(:,[n1+1:end]);301 Psit(ip)=Psit1(i);302 disp(Psit);303 end304 Psit=Psit';

305 %三、构造对应齐次方程，并解出基础解系306 b1=zeros(size(b));307

308 px=ones(1,n1);309 px(x_fqc)=0;310 px=find(px);311 npx=length(x_fqc);312 a1=a;313 %for i=1:npx314 % ipx=px(i);315 % a1(:,ipx)=a1(:,ipx)*-1;316 %disp(a1);317 %end318

319 d1=CDBH_for_sov_JZFC(a1,b1);320 %四、构造通解321 tj=[d1,Psit];322 %五、输出结果323 d=tj;324 disp '非齐次，一般解：'

325 disp 'x=k1*psi1+...+kn*psin+psit'

326 else

327 error('输入有误，无法计算');328 return;329 end330 end