gaussxqf.m

为配套王能超版本数值分析简明教程编写,方便学生教师演示!

function B=gaussxqf(A,b)%此算法的好处是可以计算矩阵的一切形式的解，并不光光是矩阵是方阵，而且行列式为零；
B=[A,b];
[n,m]=size(B);
y=[1:m-1];
q=0;
for i=1:n;
    if i>m-1;
        q=3;
        break;
    end
    a=0;
    l=i;
    while a==0;%判断是否此列i的以下的每个元素都为零，
               %若为零则a为零，此循环目的是找到有一列的元素不为零为止；
        k=i;
        p=abs(B(i,l));
        for j=i:n;%此循环目的是在此列中找到绝对值最大的元素，
                  %此好处是防止在计算机中大数吃小数的，从而导制结果与精确值相差很大
            if p<abs(B(j,l))
                p=abs(B(j,l));
                k=j;
            end
        end
        a=round(abs(B(k,l))*1000)/1000;
        if a==0;
            l=l+1;
            if l>m-1;
                q=2;
                break;
            end
        end
    end
    if q==2;
       q=1;
       break;
   end
    if a~=0;
        if k~=i;%如果k与j不相等则把矩阵B的k行与i行相互交换；
            C=eye(n);
            C(i,i)=0;C(i,k)=1;
            C(k,k)=0;C(k,i)=1;
            B=C*B;
        end
        if i~=l;%如果只有到第l列时才有元素不为零，且i与l不相等，则进行i列与l相互交换；
            D=eye(m);D1=eye(m-1);
            D(i,l)=1;D1(i,l)=1;
            D(i,i)=0;D1(i,i)=0;
            D(l,l)=0;D1(l,l)=0;
            D(l,i)=1;D1(l,i)=1;
            B=B*D;y=y*D1;
        end
        B(i,:)=B(i,:)/B(i,i);
        for o=i+1:n;%B矩阵的每一行除了第i行都用第i行的aji/aii倍去减第j行，使得每一行的第j列元素为零；
            if o~=i;
                B(o,:)=B(o,:)-B(i,:)*B(o,i);
            end
        end
    end
    B=round(B*1000000)/1000000;%此功能是为了在计算中由于计算机的计算的误差使得原来本为零的数
                               %不为零的从而导致本误解的应该跳出循环的而未能跳出
end
B=round(B*10000)/10000;%此功能是为了在计算中由于计算机的计算的误差使得原来本为零的数
                               %不为零的从而导致判断失误，本误解的而变成了有解的（值得注意的是在
                               %此当中并不精确从而导致的误差是不可避免的经过此处理后计算并非有些
                               %误差而已更加严重的是有可能导致计算的错误此算法需要改进，在此无能为力）
if q==1;%此情况是在据阵化到某一行下面的元素全部为零（其中m列出外）
    if B(i:n,m)==zeros(n-i+1,1)
        fprintf('此方程组的基础解系为：（其中自由未知量为:');
        for t=i:m-1;
            fprintf('%c%d\t','x',y(t));
        end
        fprintf('）\n');
        fprintf('解的向量顺序为：\n')
        for t=1:m-1;
            fprintf('%c%d\t','x',y(t));
        end
        fprintf('\n解的向量为:\n');
        for t=i:m-1;%根据基础解析的定义分别以t位置的元素取一其余的元素为0；
            x1(i:m-1)=0;
            x1(t)=1;
            x=zeros(1,i-1);
            for r=i-1:-1:1;%分别求每个元素在基础解析下的值
               x(r)=B(r,m)-B(r,r+1:i-1)*x(r+1:i-1)'-B(r,i:m-1)*x1(i:m-1)';
            end
            fprintf('%d\t',x);
            fprintf('%d\t',x1(i:m-1));
            fprintf('\n');
        end
    else
        disp('此方程组没有实数解1');
    end
elseif q==3;%此情况是矩阵A的行数大于列数
    if B(i:n,m)~=zeros(i:n,1);%由于此矩阵的行小于列，如果化到最后b下面的全为零则无解
        disp('此方程组没有实数解2');
    else%由于此矩阵的行小于列，如果化到最后b下面的全为零则有解
        %for r=1:i-1;
        %    x(r)=B(r,m);
        %end
        x=zeros(1,n);
        for r=i-1:-1:1;
            x(r)=B(r,m)-B(r,r+1:n)*x(r+1:n)';%先把矩阵B化为下三角矩阵，然后从上往下进行回代分别求出每一个函数的值
                                         %值得注意的是这种可能是最最常见的形式。
        end
        for t=1:i-1;
            fprintf('%c%d=%f\n','x',y(t),x((t)));
        end
    end
else%此时阵A必须是方阵，而且行列式的值不为零；
    disp('此方程组的解为：');
    fprintf('\n');
    x=zeros(1,n);
    for r=n:-1:1;
        x(r)=B(r,m)-B(r,r+1:n)*x(r+1:n)';%先把矩阵B化为下三角矩阵，然后从上往下进行回代分别求出每一个函数的值
                                         %值得注意的是这种可能是最最常见的形式。
    end
    for t=1:m-1;
        fprintf('%c%d=%f\n','x',y(t),x(t));
    end
end