基于主元素思想的householder正交法解方程组

基于主元素思想的householder正交法解方程组

• 主程序

%% 随机产生n阶实矩阵A和边值b
clc
clear
n=100
A=randi(100,n,n)
b=randi(100,n,1);
A_rec=A;%A=A_rec;
b_rec=b;%b=b_rec;
%% 正交变换
for j=1:n-1
 M=[A(j:n,:) b(j:n)];
 M=sortrows(M,-j);
 A(j:n,:)=M(:,1:n);
 b(j:n)=M(:,n+1);
 p=householder(A(j:n,j));
 unimat=eye(j-1);
 p=blkdiag(unimat,p);
 A=p*A;
 b=p*b;
end
%% 反代求值
x=zeros(n,1);
for k=n:-1:1
    x(k)=(b(k)-sum(A(k,k+1:n)*x(k+1:n)))/A(k,k);
end
%% 将计算结果于计算机计算结果对比
x
x_com=A_rec\b_rec

• 子程序householder

function p=householder(x)
l=length(x);
x=reshape(x,l,1);
sigma=norm(x,2);
abs_x1=abs(x(1));
re_x1=real(x(1));
im_x1=imag(x(1));
alpha=atan(im_x1/re_x1);
if x(1)==0 
 alpha=0;
end
if re_x1<0
if im_x1>=0
        alpha=alpha+pi;
else
        alpha=alpha-pi;
end
end
e1=[1;zeros(l-1,1)];
e_ialpha=exp(1i*alpha);(1+i)/sqrt(2);
k=-sigma*e_ialpha;
w=(x-k*e1)/sqrt(2*sigma^2+2*sigma*abs_x1);
p=eye(l)-2*(w*w');
end

householder变换也是QR分解的一个重要方法。