solvebyqr.m

来自「MATLAB科学计算与工程分析源代码源程序4」· M 代码 · 共 42 行

function x=solvebyQR(A,b)
% 该函数利用QR分解法求线性方程组Ax=b的解
% 编写日期：2007-5-14

flag=isexist(A,b); %调用第一小节中的isexist函数判断方程组解的情况
if flag==0
    disp('该方程组无解！');
    x=[];
    return;
else
    r=rank(A);
    [m,n]=size(A);
    [Q,R]=qr(A);
    b=Q'*b;
    
    % 解Rx=b得原方程组的一个特解
    x0(r)=b(r)/R(r,r);
    if r>1
        for i=r-1:-1:1
            x0(i)=(b(i)-R(i,i+1:r)*x0(i+1:r)')/R(i,i);
        end
    end
    x0=x0';
    
    if flag==1  %若方程组有唯一解
        x=x0;
        return;
    else        %若方程组有无穷多解
        format rat;
        Z=null(A,'r'); %求出对应齐次方程组的基础解系
        [mZ,nZ]=size(Z);
        x0(r+1:n)=0;
        for i=1:nZ
            t=sym(char([107 48+i]));
            k(i)=t;      %取k=[k1,...,kr];            
        end
        x=x0;         
        for i=1:nZ          
            x=x+k(i)*Z(:,i); %将方程组的通解表示为特解加对应齐次通解形式
        end        
    end
end