线性方程组的迭代解法-豆柴文库

您所在位置：网站首页 / 线性方程组的迭代解法.doc / 文档详情

在线预览结束，喜欢就下载吧，查找使用更方便

10 金币

下载文档

/ 3

下载提示

如果您无法下载资料，请参考说明：

1、部分资料下载需要金币，请确保您的账户上有足够的金币

2、已购买过的文档，再次下载不重复扣费

3、资料包下载后请先用软件解压，在使用对应软件打开

文本预览

实验三线性方程组的迭代求解法
题目：求解方程组AX=b其中A=[-4111;1-411;11-41;111-4];b=[1111]';
方法：Jacobi迭代法，Gauss-Seidel迭代法
程序：
（1）Jacobi迭代法
function[x,k,flag]=jacobi(A,b,delta,max1)
ifnargin<4max1=100;end
ifnargin<3delta=1e-4;end
n=length(A);k=0;
x=zeros(n,1);y=zeros(n,1);flag='ok';
while1
fori=1:n
y(i)=b(i);
forj=1:n
ifj~=i
y(i)=y(i)-A(i,j)*x(j);
end
end
ifabs(A(i,i))<1e-10|k==max1
flag='fail';return;
end
y(i)=y(i)/A(i,i);
end
ifnorm(y-x,inf)<delta
break;
end
x=y;k=k+1;
end
（2）Gauss-Seidel迭代法
function[x,k,flag]=gao(A,b,delta,max1)
ifnargin<4max1=100;end
ifnargin<3delta=1e-4;end
n=length(A);k=0;
x=zeros(n,1);y=zeros(n,1);flag='ok';
while1
y=x;
fori=1:n
z=b(i);
forj=1:n
ifj~=i
z=z-A(i,j)*x(j);
end
end
ifabs(A(i,i))<1e-10|k==max1
flag='fail';return;
end
z=z/A(i,i);
x(i)=z;
end
ifnorm(y-x,inf)<delta
break;
end
k=k+1;
end
实验结果：
（1）Jacobi迭代法
>>A=[-4111;1-411;11-41;111-4];
>>b=[1111]';
>>formatcompact
>>[x,k,flag]=jacobi(A,b)
x=
-0.99968252072859
-0.99968252072859
-0.99968252072859
-0.99968252072859
k=
28
flag=
ok
（2）Gauss-Seidel迭代法
>>A=[-4111;1-411;11-41;111-4];
>>b=[1111]';
>>formatlong
>>[x,k,flag]=gao(A,b)
x=
-0.99990189431986
-0.99991475824547
-0.99992593541266
-0.99993564699450
k=
16
flag=
ok
五、拓展：希尔伯特矩阵，这里只取n=4时的情况。
>>H4=[11/21/31/4;1/21/31/41/5;1/31/41/51/6;1/41/51/61/7]
H4=
1.00000.50000.33330.2500
0.50000.33330.25000.2000
0.33330.25000.20000.1667
0.25000.20000.16670.1429
>>cond(H4,inf)
ans=
2.8375e+004
>>b=[2.08331.28330.95000.7595]'
b=
2.0833
1.2833
0.9500
0.7595
>>A=H4
A=
1.00000.50000.33330.2500
0.50000.33330.25000.2000
0.33330.25000.20000.1667
0.25000.20000.16670.1429
>>[x,k,flag]=jacobi(A,b)
x=
1.0e+041*
-0.8411
-1.6834
-2.1499
-2.4542
k=
100
flag=
fail
>>[x,k,flag]=gao(A,b)
x=
1.0018
0.9564
1.1291
0.9070
k=
100
flag=
fail
由以上结果可知，雅可比迭代法、高斯-塞德尔迭代法都不能求解出方程组。而由的无穷范数条件数为28375可知是一个病态矩阵，这个方程组本身是个病态问题，和算法没有关系。