c3_4.m

这是我从一本基于MATLAB实现的数值分析算法的书中自带的光盘程序

%C3_4.m
%不定方程组的求解
%A：系数矩阵；d：右端常数项；a0：增广矩阵；a:选主元后的增广矩阵

function C3_4
A=[2 3 1 4 1;2 3 1 1 -1;4 6 -1 1 2];
d=[6 1 5]';
a0=[A,d];
m=length(a0(:,1));
n=length(a0(1,:));
g=length(a0(:,1));
%选主元
P=eye(g);
for k=1:g 
    for p=k:g
        ekp=zeros(g);
        I=eye(g);
        ekp=I-(I(:,k)-I(:,p))*((I(:,k)-I(:,p))');
       if abs(a0(p,k))==max(abs(a0(k:g,k)))
          P=ekp*P;
       end
    end
end
P;
a=P*a0;
%用约当消去法
    a(1,:)=a(1,:)/a(1,1)
    for j=2:m
        a(j,:)=a(j,:)+a(1,:)*(-a(j,1));
    end
    for i=2:m
        for j=i:m+1
            if a(i,j-1)==0 
                if a(i,j)~=0
                   a(i,:)=a(i,:)/a(i,j);
                   for l=1:i-1
                       a(l,:)=a(l,:)+(-a(l,j)*a(i,:));
                   end
                   for l=i+1:m
                       a(l,:)=a(l,:)+(-a(l,j)*a(i,:));
                   end
                end
            end
        end
    end
    fprintf('最后的增广矩阵为；\n')
    a
fprintf('u= %6.5f + (%6.5f) v + (%6.5f) y \n',a(1,n),-a(1,2),-a(1,n-1));
fprintf('w= %6.5f + (%6.5f) y \n',a(2,n),-a(2,n-1));
fprintf('x= %6.5f + (%6.5f) y \n',a(3,n),-a(3,n-1));
%t=input('自由变量【v y】：');