选主元LU分解
实验内容:列选主元LU分解和全选主元LU分解求解线性方程组
计算方法:
- 全选主元消元法
1.1 初始化
根据参数A、b,记录下矩阵、右端项的尺寸n;
以得到的尺寸n初始化解向量x;
同时,以尺寸(n-1)初始化一个记录列交换顺序的向量tcolum;
1.2 进入循环
for k=1:n-1
选主元:|A(p1,p2)|=max{|A(i,j)|:i=k:n,j=k:n}
记录列交换顺序:tcolum(k)=p2
换行:A(k,k:n)与A(p1,k:n)互换,同时b(k)与b(p1)互换
换列:A(k:n,k)与A(k:n,p2)互换
消元:
~~~ for i=k+1:n
~~~~~~ A(i,k)=A(i,k)/A(k,k);
~~~~~~ for j=k+1:n
~~~~~~~~~ A(i,j)=A(i,j)-A(i,k)A(k,j)
~~~~~~ end
~~~ end
end
注:记录列交换顺序,举例:tcolum(2)=5表示第二列和第五列互换。
寻找第k个主元时,第k个主元原始位置的列数总要大于k。
最多需要交换n-1次(矩阵A维数nn),所以向量tcolum是n-1维的。
1.3 回代求解
x(n)=b(n)/A(n,n)
for k=n-1: -1:1
~~~~ s=sum(A(k,k+1:n));
~~~~ x(k)=(b(k)-s)/A(k,k);
end
1.4 理顺解向量
for i=n-1: -1:1
~~~~ 交换x(i)与x(tcolum(i))
end - 列选主元消元法
2.1 初始化
根据参数A、b,记录下矩阵、右端项的尺寸n;
以得到的尺寸n初始化解向量x;
2.2 进入循环
for k=1:n-1
~~~~ 选列主元:|A(p1,k)|=max{|A(i,j)|:i=k:n}
~~~~ 换行:A(k,k:n)与A(p1,k:n)互换,同时b(k)与b(p1)互换
~~~~ 消元:
~~~~ for i=k+1:n
~~~~~~~~ A(i,k)=A(i,k)/A(k,k);
~~~~~~~~ for j=k+1:n
~~~~~~~~~~~~ A(i,j)=A(i,j)-A(i,k)*A(k,j)
~~~~~~~~ end
~~~~ end
end
2.3 回代求解
x(n)=b(n)/A(n,n)
for k=n-1: -1:1
~~~~ s=sum(A(k,k+1:n));
~~~~ x(k)=(b(k)-s)/A(k,k);
end
计算结果:
全选主元计算结果的残量范数为8.4428e-011
列选主元计算结果的残量范数为1.8175e-010
源程序:
1.定义全选主元函数
function x=Lu_Decomplete(A,b)
%定义全选主元--解Ax=b
[n,n2]=size(A);
n3=length(b);
if n~=n2|n~=n3 %判断输入的合法性
error('wrong input!');
end
x=ones(n,1);%记录解
tcolum=zeros(n-1,1);%记录列变换
for k=1:n-1 %从(1,1)位置到(n-1,n-1)位置
%进行全选主元
p1=k;p2=k;p_max=abs(A(k,k));%p1、p2分别记录/将选主元/的行、列
for i=k:n
for j=k:n
if p_max<abs(A(i,j))
p_max=abs(A(i,j));
p1=i;p2=j;
end
end
end
%对A、b换行
if p1>k
for i=k:n
t=A(k,i);
A(k,i)=A(p1,i);
A(p1,i)=t;
end
t=b(k);
b(k)=b(p1);
b(p1)=t;
end
tcolum(k)=p2;
%对A换列
if p2>k
for i=1:n
t=A(i,k);
A(i,k)=A(i,p2);
A(i,p2)=t;
end
end
%消元
if A(k,k)<eps
error('Fail solving');
end
for i=k+1:n
A(i,k)=A(i,k)/A(k,k);
for j=k+1:n
A(i,j)=A(i,j)-A(i,k)*A(k,j);
end
b(i)=b(i)-A(i,k)*b(k);
end
end
%回代求解
x(n)=b(n)/A(n,n);
for i=n-1:-1:1
z=0;
for j=i+1:n
z=z+x(j)*A(i,j);
end
x(i)=(b(i)-z)/A(i,i);
end
%理顺解向量!!!重点在此!!
for i=n-1:-1:1
t=x(i);
x(i)=x(tcolum(i));
x(tcolum(i))=t;
end
2.定义列选主元函数
function x=Lu_Decolum(A,b)
%定义列选主元LU分解
[n,n2]=size(A);
n3=length(b);
if n~=n2|n~=n3 %判断输入的合法性
error('wrong input!');
end
x=ones(n,1);%初始化解向量
for k=1:n-1
%选出列主元
p=k;p_max=abs(A(k,k));
for i=k+1:n
if p_max<abs(A(i,k))
p_max=abs(A(i,k));
p=i;
end
end
%据选出的列主元进行换行
if p>k
for i=k:n
t=A(k,i);
A(k,i)=A(p,i);
A(p,i)=t;
end
t1=b(k);
b(k)=b(p);
b(p)=t1;
end
%对矩阵进行LU分解
if A(k,k)==0
break;
end
for j=k+1:n
m=A(j,k)/A(k,k);
for i=k+1:n
A(j,i)=A(j,i)-m*A(k,i);
end
b(j)=b(j)-m*b(k);
end
end
%回代法解方程组
x(n)=b(n)/A(n,n);
for k=n-1:-1:1
z=0;
for i=k+1:n
z=z+x(i)*A(k,i);
end
x(k)=(b(k)-z)/A(k,k);
end
3.调用函数计算结果
在命令行中,
输入A=randn(1000,1000);b=ones(1000,1);
调用全选主元函数:
x1=Lu_Decomplete(A,b);
计算此时结果的残量范数:
norm(Ax1-b,2).
同理,调用列选主元函数:
x2=Lu_Decolum(A,b);
计算此时结果的残量范数:
norm(Ax2-b,2).