数值分析Hilbert矩阵病态线性方程组的求解Matlab程序.docx
- 文档编号:15512625
- 上传时间:2023-07-05
- 格式:DOCX
- 页数:12
- 大小:359.43KB
数值分析Hilbert矩阵病态线性方程组的求解Matlab程序.docx
《数值分析Hilbert矩阵病态线性方程组的求解Matlab程序.docx》由会员分享,可在线阅读,更多相关《数值分析Hilbert矩阵病态线性方程组的求解Matlab程序.docx(12页珍藏版)》请在冰点文库上搜索。
数值分析Hilbert矩阵病态线性方程组的求解Matlab程序
(Hilbert矩阵)病态线性方程组的求解
理论分析表明,数值求解病态线性方程组很困难。
考虑求解如下的线性方程组的求解Hx=b,期中H是Hilbert矩阵,
,
,i,j=1,2,…,n
1.估计矩阵的2条件数和阶数的关系
2.对不同的n,取
,分别用Gauss消去,Jacobi迭代,Gauss-seidel迭代,SOR迭代和共轭梯度法求解,比较结果。
3.结合计算结果,试讨论病态线性方程组的求解。
第1小题:
condition.m%第1小题程序
t1=20;%阶数n=20
x1=1:
t1;
y1=1:
t1;
fori=1:
t1
H=hilb(i);
y1(i)=log(cond(H));
end
plot(x1,y1);
xlabel('阶数n');
ylabel('2-条件数的对数(log(cond(H))');
title('2-条件数的对数(log(cond(H))与阶数n的关系图');
t2=200;%阶数n=200
x2=1:
t2;
y2=1:
t2;
fori=1:
t2
H=hilb(i);
y2(i)=log(cond(H));
end
plot(x2,y2);
xlabel('阶数n');
ylabel('2-条件数的对数(log(cond(H))');
title('2-条件数的对数(log(cond(H))与阶数n的关系图');
画出Hilbert矩阵2-条件数的对数和阶数的关系
n=200时
n=20时
从图中可以看出,
1)在n小于等于13之前,图像近似直线
log(cond(H))~1.519n-1.833
2)在n大于13之后,图像趋于平缓,并在一定范围内上下波动,同时随着n的增加稍有上升的趋势
第2小题:
solve.m%m第2小题主程序
N=4000;
xGauss=zeros(N,1);
xJacobi=zeros(N,1);
xnJ=zeros(N,1);
xGS=zeros(N,1);
xnGS=zeros(N,1);
xSOR=zeros(N,1);
xnSOR=zeros(N,1);
xCG=zeros(N,1);
xnCG=zeros(N,1);
forn=1:
N;
x=ones(n,1);
t=1.1;%初始值偏差
x0=t*x;%迭代初始值
e=1.0e-8;%给定的误差
A=hilb(n);
b=A*x;
max=100000000000;%可能最大的迭代次数
w=0.5;%SOR迭代的松弛因子
G=Gauss(A,b);
[J,nJ]=Jacobi(A,b,x0,e,max);
[GS,nGS]=G_S(A,b,x0,e,max);
[S_R,nS_R]=SOR(A,b,x0,e,max,w);
[C_G,nC_G]=CG(A,b,x0,e,max);
normG=norm(G'-x);
xGauss(n)=normG;
normJ=norm(J-x);
nJ;
xJacobi(n)=normJ;
xnJ(n)=nJ;
normGS=norm(GS-x);
nGS;
xGS(n)=normGS;
xnGS(n)=nGS;
normS_R=norm(S_R-x);
nS_R;
xSOR(n)=normS_R;
xnSOR(n)=nS_R;
normC_G=norm(C_G-x);
nC_G;
xCG(n)=normC_G;
xnCG(n)=nC_G;
end
Gauss.m%Gauss消去法
functionx=Gauss(A,b)
n=length(b);
l=zeros(n,n);x=zeros(1,n);
%消去过程
fori=1:
n-1
forj=i+1:
n
l(j,i)=A(j,i)/A(i,i);
fork=i:
n
A(j,k)=A(j,k)-l(j,i)*A(i,k);
end
b(j)=b(j)-l(j,i)*b(i);
end
end
%回代过程
x(n)=b(n)/A(n,n);
fori=n-1:
-1:
1
c=A(i,:
).*x;
x(i)=(b(i)-sum(c(i+1:
n)))/A(i,i);
end
Jacobi.m%Jacobi迭代,x0表示迭代初值,e表示允许误差(迭代停止条件),n表示迭代次数,m可能最大的迭代次数
function[x,n]=Jacobi(A,b,x0,e,m)
n=length(A);
D=diag(diag(A));
U=-triu(A,1);
L=-tril(A,-1);
B=D\(L+U);
f=D\b;
x=B*x0+f;
n=1;
whilenorm(x-x0)>e
x0=x;
x=B*x0+f;
n=n+1;
ifn>m
disp('Jacobi迭代次数过多,迭代可能不收敛');
break;
end
end
G_S.m%Gauss-Seidel迭代,x0表示迭代初值,e表示允许误差(迭代停止条件),n表示迭代次数,m可能最大的迭代次数
function[x,n]=G_S(A,b,x0,e,m)
n=length(A);
D=diag(diag(A));
U=-triu(A,1);
L=-tril(A,-1);
B=(D-L)\U;
f=(D-L)\b;
x=B*x0+f;
n=1;
whilenorm(x-x0)>e
x0=x;
x=B*x0+f;
n=n+1;
ifn>m
disp('Gauss-Seidel迭代次数过多,迭代可能不收敛');
break;
end
end
SOR.m%SOR超松弛迭代,x0表示迭代初值,e表示允许误差(迭代停止条件),n表示迭代次数,m可能最大的迭代次数,w松弛因子
function[x,n]=SOR(A,b,x0,e,m,w)
n=length(A);
D=diag(diag(A));
U=-triu(A,1);
L=-tril(A,-1);
B=(D-w*L)\((1-w)*D+w*U);
f=(D-w*L)\b*w;
x=B*x0+f;
n=1;
whilenorm(x-x0)>e
x0=x;
x=B*x0+f;
n=n+1;
ifn>m
disp('SOR超松弛迭代次数过多,迭代可能不收敛');
break;
end
end
CG.m%CG共轭梯度法,x0表示迭代初值,e表示允许误差(迭代停止条件),n表示迭代次数,m可能最大的迭代次数
function[x,n]=CG(A,b,x0,e,m)
r=b-A*x0;
p=r;
alpha=(r'*r)/(p'*(A*p));
x=x0+alpha*p;
r1=b-A*x;
n=1;
whilenorm(r1)>e
belta=(r1'*r1)/(r'*r);
p=r1+belta*p;
r=r1;
x0=x;
alpha=(r'*r)/(p'*(A*p));
x=x0+alpha*p;
r1=b-A*x;
n=n+1;
ifn>m
disp('CG共轭梯度法迭代次数过多,迭代可能不收敛');
break;
end
end
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数值 分析 Hilbert 矩阵 病态 线性方程组 求解 Matlab 程序