Hilbert代数方程组的数值解法.docx
- 文档编号:11517597
- 上传时间:2023-03-02
- 格式:DOCX
- 页数:16
- 大小:43.13KB
Hilbert代数方程组的数值解法.docx
《Hilbert代数方程组的数值解法.docx》由会员分享,可在线阅读,更多相关《Hilbert代数方程组的数值解法.docx(16页珍藏版)》请在冰豆网上搜索。
Hilbert代数方程组的数值解法
Hilbert代数方程组的数值解法
1Hilbert矩阵的条件数和矩阵的阶数的关系
编制Matlab程序:
clear,clc
formatlong
%formatshort
forn=1:
14;
A=hilb(n);
%A=rand(n);
condA(n)=cond(A,inf);
end
plot(log10(condA))
得出图形
图1:
Hilbert矩阵和阶次的关系
由图可见a.Hilbert矩阵的2-条件数的对数与阶次近似正比;
b.阶次超过12后,计算困难,舍入误差会导致计算不准
结论:
Hilbert矩阵的2-条件数随阶次指数增长,关系大概是:
Log10(Cond(A))=1.33791720780254(n-1)
2.各种求解方法的对比
2.1直接法:
Gauss消去方法(程序见附录,下同)
在双精度型变量下,即eps=2.220446049250313e-016,对于直接法Gauss消去方法,求解Hilbert矩阵方程组,在阶次n=13时,解得:
X=[10.999971.00130.978811.1906-0.0354414.6171-7.396714.089-12.5429.9162-2.38171.5624]
出现明显错误,可见Gauss消去方法对病态问题比较敏感。
求解阶次不高。
2.2Jacobi迭代法
很遗憾,jakobi迭代法的迭代矩阵的谱半径随阶次线性上升(程序见附录2),普遍大于1,计算结果发送,无法迭代出满意的结果。
需放弃
图2:
Jacobi迭代矩阵的谱半径
2.3Gauss-Seidel迭代与SOR迭代求解
先分析SOR迭代矩阵的谱半径和松弛因子w的关系
图3:
SOR迭代矩阵的谱半径和松弛因子w的关系
可以发现SOR迭代和G-S迭代的谱半径都普遍接近于1,因此松弛因子的选择对计算结果的影响不大。
我们选择w=1,即G-S迭代来做数值实验:
选定e=0.000001;
迭代次数:
15610次
下降曲线:
图4:
G-S迭代下降曲线
取w=1.5时,
迭代次数需58798次
前100步下降曲线
图5:
SOR(w=1.5)下降曲线
w=0.5迭代次数16290
图6:
SOR(w=0.5)下降曲线
2.4SOR按照寻优算法的迭代求解
寻优得出w=0.1时,谱半径相对较小,迭代次数只需3689次,下降曲线也非常快,如图7
图7:
SOR(w=0.1)下降曲线
2.5共轭梯度法迭代求解
奇迹发生了,迭代只需要9次就完成了,下降也非常快,见图
图8:
共轭梯度法下降曲线
2.6各算法综合对比
综上所述,在双精度的数值精度下和e=0.00001下,在笔记本可接受的计算时间里,做各种算法的对比如下:
Gauss
Jacobi
G-S
SOR
CG
能计算到n阶
12
2
1200
2000
5500
迭代次数
/
/
最多
中等
较少
计算时间
/
慢
较快
超快
稳定性
较差
最差
一般
较好
很好
3.讨论病态方程组的求解
3.1对于Hilbert病态方程组的求解
Gauss直接法和Jacobi迭代法都比较差,G-S法和SOR虽能求解,但由于其条件数都接近于1,收敛太慢,计算时间太长,因此实际可计算的阶次不高。
3.2采用高精度计算
我们分别采用单精度和双精度型变量,以Gauss消去方法为求解办法,做数值实验,
x2=gauss(A,b)
x3=gauss(single(A),single(b))
当阶次n=6时,得
x2=
0.99999999999908
1.00000000002630
0.99999999982172
1.00000000046435
0.99999999948691
1.00000000020235
x3=
1.0003130
0.9906988
1.0645127
0.8293791
1.1905806
0.9242350
可以发现,x3已经明显偏离了正解,得出结论,采用高精度计算所得解有更多有效位。
3.3采用预处理方法
预处理共轭梯度法见附录,采用后,可以发现,计算4000阶的问题,速度加快很多,迭代次数也减少了。
附录1Gauss消去方法
functionx=Gauss(a,b)
%Gauss消去法
[n,m]=size(a);
nb=length(b);
det=1;%存储行列式值
x=zeros(n,1);
fork=1:
n-1
fori=k+1:
n
ifa(k,k)==0
return
end
m=a(i,k)/a(k,k);
forj=k+1:
n
a(i,j)=a(i,j)-m*a(k,j);
end
b(i)=b(i)-m*b(k);
end
det=det*a(k,k);
end
det=det*a(n,n);
fork=n:
-1:
1%回代
forj=k+1:
n
b(k)=b(k)-a(k,j)*x(j);
end
x(k)=b(k)/a(k,k);
end
附录2Jacobi迭代的收敛性分析
clear,clc
formatlong
%formatshort
forn=1:
25;
A=hilb(n);
D=diag(diag(A));
B=inv(D)*(D-A);
x=ones(n,1);
e=0.000001;
e=eps;
b=A*x;
f=inv(D)*b;
%%谱半径大于等于1就不收敛
p(n)=max(abs(eig(B)));
end
plot(p)
xlabel('阶次')
ylabel('谱半径')
title('Jacobi迭代收敛性分析')
附录3G-S迭代
function[x,k,Ve]=gauss_seidel(A,b,e,M)
[NAr,NA]=size(A);
ifNAr~=NA,error('Aisnotasquarematrix');end
%检查A,b的大小是否匹配
[Nbr,Nb]=size(b);
if((Nbr~=NA)&&(1~=Nb)),error('Aandbhavenon-compatibledimensions');end
%e不能小于计算机的最小精度
ife x=zeros(size(b));%初始解设置为与b同型的零向量 k=0;%迭代次数的记数变量,初始量设为0 r=1+e;%前后项之间的误差 %构造D、B、f矩阵 D=diag(diag(A)); U=-triu(A,1); L=-tril(A,-1); DsubL=(D-L); B=inv(DsubL)*U; f=inv(DsubL)*b; %为了画图,返回误差向量 Ve=zeros(100,1); %谱半径大于等于1就不收敛 p=max(abs(eig(B))); ifp>=1 sprintf('G-Siterationmethodnotconvergent') %return end while(r>=e)&&(k~=M) x0=x; x=B*x0+f; k=k+1; r=max(abs(x-x0)); % %ifk<=size(Ve,1) %Ve(k,1)=r; %end end 附录4SOR迭代 function[x,n]=SOR(A,b,w,e,M) x0=zeros(size(b));%初始解设置为与b同型的零向量 %if(w<=0||w>=2) %error; %return; %end D=diag(diag(A));%求A的对角矩阵 L=-tril(A,-1);%求A的下三角阵 U=-triu(A,1);%求A的上三角阵 B=inv(D-L*w)*((1-w)*D+w*U); f=w*inv((D-L*w))*b; n=0;%迭代次数 r=1+e; whiler>=e x=B*x0+f; n=n+1; r=max(abs(x-x0)); %Ve(n)=r; x0=x; if(n>=M) disp('Warning: 迭代次数太多,可能不收敛! '); return; end end %figure %stem(Ve(1: 100)) %title('SOR迭代误差曲线'); %xlabel('迭代次数'); %ylabel('误差');... 附录5CG方法 function[x,n]=conjgrad(A,b,x0) r1=b-A*x0; p=r1; n=0; fori=1: rank(A)%以下过程可参考算法流程 if(dot(p,A*p)<1.0e-50)%循环结束条件 break; end alpha=dot(r1,r1)/dot(p,A*p); x=x0+alpha*p; r2=r1-alpha*A*p; if(r2<1.0e-50)%循环结束条件 break; end belta=dot(r2,r2)/dot(r1,r1); p=r2+belta*p; n=n+1; end 附录6PCG方法 function[x,n]=preconjgrad(A,b,M,eps) x0=zeros(size(b));%初始解设置为与b同型的零向量 ifnargin==4 eps=1.0e-6; end r1=b-A*x0; iM=inv(M); z1=iM*r1; p=z1; n=0; tol=1; whiletol>=eps alpha=dot(r1,z1)/dot(p,A*p); x=x0+alpha*p; r2=r1-alpha*A*p; z2=iM*r2; belta=dot(r2,z2)/dot(r1,z1); p=z2+belta*p; n=n+1; tol=norm(x-x0); x0=x;%更新迭代值 r1=r2; z1=z2; end
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- Hilbert 代数 方程组 数值 解法