Gauss列主元消去法、QR(MATLAB)

例：用Gauss列主元消去法、QR方法求解如下方程组：

22 41

4 2 2312 x1 1 3 1 x2 2 .01 x3 1 23 x4 0

1. 1）Gauss列主元法源程序：

function x=Gauss(A,b)

[m,n]=size(A);

if m~=n

error('矩阵不是方阵')

return

end

B=[A,b];

n=length(A);

for j=1:n-1

q=[zeros(j-1,1);B(j:n,j)];

[c,r]=max(abs(q)); %c为列主元，r为所在行 if r~=j

temp=B(j,:); %交换两行

B(j,:)=B(r,:);

B(r,:)=temp;

end

for i=j+1:n

B(i,:)=B(i,:)-B(j,:)*(B(i,j)/c);

end

end

x(n)=B(n,n+1)/B(n,n);

for i=n-1:-1:1

for j=i:n-1

B(i,n+1)=B(i,n+1)-B(i,j+1)*x(j+1);

end

x(i)=B(i,n+1)/B(i,i);

end

2)在命令窗口输入A，b,得到x的近似解：

>> A=[2,2,1,2;4,1,3,-1;-4,-2,0,1;2,3,2,3];

>> b=[1;2;1;0];

>> x=Gauss(A,b)

1.5417

-2.7500

0.0833

1.6667

2. 1）QR方法源程序：

function [Q,R,X]=qrfj(A,b)

[m,n]=size(A);

if m<n

error('A不符合规则')

return

end

R=A;

Q=eye(n);

for i=1:n-1

a=R(i:n,i);

e=[1;zeros(n-i,1)];

w=a-norm(a)*e;

Hw=eye(n-i+1)-(2/(w'*w))*(w*w');

H=blkdiag(eye(i-1),Hw); Q=Q*H;

R=H*R;

end

Y=Q'*b;

X(n)=Y(n)/R(n,n);

for i=n-1:-1:1

for j=i:n-1

Y(i)=Y(i)-R(i,j+1)*X(j+1);

end

X(i)=Y(i)/R(i,i);

end

2)在命令窗口输入A，b,得到x的近似解：

>> A=[2,2,1,2;4,1,3,-1;-4,-2,0,1;2,3,2,3];

>> b=[1;2;1;0];

>> [Q,R,X]=qrfj(A,b)

Q =

0.3162 0.3705 -0.0284 -0.8729

0.6325 -0.4940 0.5966 -0.0000

-0.6325 0.0823 0.7386 -0.2182

0.3162 0.7822 0.3125 0.4364 块对角矩阵 %

6.3246 3.4785 2.8460 0.3162 0.0000 2.4290 0.4529 3.6641 0.0000 0.0000 2.3864 1.0227

-0.0000 -0.0000 -0.0000 -0.6547

X =

1.5417 -2.7500 0.0833 1.6667