MATLAB中矩陣方程求解的實現
一、矩陣方程
1、定義：
2、分類
http://naotu.baidu.com/file/14d36860667d356a54490320cdab2950?token=56a299fe0e0815fe
二、M代碼實現
1、M代碼
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~=m2 21 error('輸入有誤，無法計算'); 22 return; 23 end 24 switch pj 25 %這種情況為其次方程組 26 case 0 27 switch rc 28 %這種情況只有零解 29 case n1 30 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行及以下第一個非零元素的位置k 63 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 end 77 end 78 end 79 i0=i1+i; 80 x_qc=[x_qc,i0]; %將此時的列位置進行記錄 81 %Step2 檢查階梯元素是否為1，若不為1，則將其化為1 82 if c(i,i0)~=0&&c(i,i0)~=0; 83 c(i,:)=c(i,:)/c(i,i0); 84 end 85 %Step3 將i0列，第i行一下的元素消減為0 86 for j=i+1:m1 87 if c(j,i0)~=0 88 c(j,:)=c(j,:)-c(j,i0)*c(i,:); 89 end 90 end 91 %Step4 將i0列，第i行一上的元素消減為0 92 for j=1:i-1 93 if c(j,i0)~=0 94 c(j,:)=c(j,:)-c(j,i0)*c(i,:); 95 end 96 end 97 end 98 %更新！檢查最后一行 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 end 105 end 106 end 107 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; %其中線性無關向量位置設為0 115 no_x_qc=find(no_x_qc); %找到非零元素位置，命名為no_x_qc 116 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 end 127 %Step5 將基礎解系中，當前列對應的位置的元素賦值為1 128 jcjx(no_x_qc(i),i)=1; 129 end 130 %disp(jcjx) 131 %三、輸出結果 132 disp '齊次型，結果為基礎解系' 133 d=jcjx; 134 disp 'x1' 135 disp '... = k1*psi1+...+kr*psir' 136 disp 'xn' 137 otherwise 138 error('輸入有誤，無法計算'); 139 return; 140 end 141 %這種話情況為非齊次方程 142 otherwise 143 %這種情況無解 144 if ra<rc 145 error('非齊次方程組無解'); 146 return; 147 %這種情況唯一解 148 elseif ra==rc&&rc==n1 %這時，a必然為方陣，但是可以使用上述齊次方程的解法 149 %思路： 150 %一、矩陣變換 151 %二、求解 152 %三、輸出結果 153 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行及以下第一個非零元素的位置k 172 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 end 186 end 187 end 188 i0=i1+i; 189 x_qc=[x_qc,i0]; %將此時的列位置進行記錄 190 %Step2 檢查階梯元素是否為1，若不為1，則將其化為1 191 if c(i,i0)~=0&&c(i,i0)~=0; 192 c(i,:)=c(i,:)/c(i,i0); 193 end 194 %Step3 將i0列，第i行一下的元素消減為0 195 for j=i+1:m1 196 if c(j,i0)~=0 197 c(j,:)=c(j,:)-c(j,i0)*c(i,:); 198 end 199 end 200 %Step4 將i0列，第i行一上的元素消減為0 201 for j=1:i-1 202 if c(j,i0)~=0 203 c(j,:)=c(j,:)-c(j,i0)*c(i,:); 204 end 205 end 206 end 207 %更新！檢查最后一行 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 end 214 end 215 end 216 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<n1 228 %思路： 229 %一、矩陣變換（利用齊次方程的代碼） 230 for i=1:m1-1; %需要對行數-1行進行行變換 231 %選中了第i行 232 %如果都為非線性相關的向量，則階梯的行列數相等，若存在線性相 233 %的向量，則需要構造一個i0來表示階梯對應的列，並用i1表示列於 234 %行數的差值。 235 i0=i1+i; 236 %因為i1的值根據本行元素具體情況確定，因此需要注意，應當在這 237 %個for循環內完成下三角、化為1、上三角的所有操作。 238 %Step01 檢查階梯元素是否等於零 239 % 若等於零，則需要與第一個不為零的行對調，若全為零， 240 % 則i1+1，並再次循環。 241 pj_01=0;%初始化以下循環的判據 242 while pj_01==0 %利用判據pj_01來判斷是否需要再次循環，在循環中，通過修改pj的值來跳出循環。 243 pj_01=1; %沒有特殊情況執行完跳出循環 244 i0=i+i1; 245 if c(i,i0)==0 246 %01 找到第i0列，i行及以下第一個非零元素的位置k 247 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 end 261 end 262 end 263 i0=i1+i; 264 x_qc=[x_qc,i0]; %將此時的列位置進行記錄 265 %Step2 檢查階梯元素是否為1，若不為1，則將其化為1 266 if c(i,i0)~=0&&c(i,i0)~=0; 267 c(i,:)=c(i,:)/c(i,i0); 268 end 269 %Step3 將i0列，第i行一下的元素消減為0 270 for j=i+1:m1 271 if c(j,i0)~=0 272 c(j,:)=c(j,:)-c(j,i0)*c(i,:); 273 end 274 end 275 %Step4 將i0列，第i行一上的元素消減為0 276 for j=1:i-1 277 if c(j,i0)~=0 278 c(j,:)=c(j,:)-c(j,i0)*c(i,:); 279 end 280 end 281 end 282 %更新！檢查最后一行 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 end 289 end 290 end 291 292 %成功化為標准型 293 %disp(c) 294 %二、求特解Psit 295 x_fqc=x_qc; %為了便於閱讀，使用x_fqc代替x_qc 296 Psit=zeros(1,n1); %初始化特解（長度：a的列數） 297 nx=length(x_fqc); %nx為線性無關向量的數量 298 for i=1:nx 299 ip=x_fqc(i); 300 Psit1=c(:,[n1+1:end]); 301 Psit(ip)=Psit1(i); 302 disp(Psit); 303 end 304 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:npx 314 % ipx=px(i); 315 % a1(:,ipx)=a1(:,ipx)*-1; 316 % disp(a1); 317 %end 318 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 end 330 end
2、例子（用於驗證）
（1）齊次方程（無限解）：
1 a=[1,1,-1,-1;2,-5,3,2;7,-7,3,1]; 2 b=[0,0,0]';
結果：
（2）非齊次方程（唯一解）：
1 a=[89,14,81,19;95,25,24,25;54,84,92,61;13,25,34,47]; 2 b=[436,317,742,353]';
結果：
>>x =
1
2
3
4
（3）非齊次方程（無限解）：
1 a=[1,-1,-1,1;1,-1,1,-3;1,-1,-2,3]; 2 b=[0,1,-1/2]';
結果：
