计算方法作业第八章

1. 考虑线性方程组Ax b =，其中系数矩阵和右端向量分别为

2111210,121

..................0121n n n

A b ⨯-⎡⎤⎡⎤⎢⎥⎢⎥--⎢⎥⎢⎥

⎢⎥⎢⎥==--⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥-⎣⎦⎣⎦ 分别采用Jacobi 迭代，Gauss-Seidel 迭代和SOR 迭代法求解，初始向量均取为零向量。 （1） 对三种方法分别进行不同停机标准时迭代次数的比较【绝对误差，相对误差】【残量，

相邻差量，下界估计方式】，误差指标为6

10Tol -=，可采用2范数或无穷范数；

（2） 绘制三种迭代方法的误差下降曲线以及残量的下降曲线（采用对数坐标系）； （3） 观测SOR 迭代方法中ω的不同取值对迭代次数的影响，并数值确定最佳opt

ω的值；

（4） 研究三种迭代算法（取最佳因子）的迭代次数与矩阵阶数倒数的关系； （5） 数值估算迭代矩阵的谱半径；

（1） J acobi ：

format long

n=input('请输入n 的值：'); A=zeros(n); for i=1:n-1 A(i,i)=2; A(i+1,i)=-1; A(i,i+1)=-1; end A(n,n)=2; b=[1]; for i=1:n-2 b=[b;0]; end b=[b;1];

x0=(linspace(0,0,n))'; tol=10^(-6); A; b; a=[]; for k=1:10000

x=[];

for i=1:n

s=0;

for j=1:n

if j~=i

s=s+A(i,j)*x0(j);

end

end

x(i)=(b(i)-s)/A(i,i);

end

aa=0;

for i=1:n

aa=aa+((x(i)-x0(i))^2);

end

a(k)=sqrt(aa);

if a(k)<=tol%判定系数为a(k)

x

k

break

else

x0=x;

end

end

K=linspace(1,k,k);

figure

semilogy(K,a,'-');

ylabel('误差')

xlabel('迭代次数')

Gauss-Seidel：

format long

n=input('请输入n的值：');

p=input('请选择判定标准代号：(1:相邻差量绝对误差；2：相邻差量相对误差；3：残量绝对误差；4：残量相对误差)');

A=zeros(n);

for i=1:n-1

A(i,i)=2;

A(i+1,i)=-1;

A(i,i+1)=-1;

end

A(n,n)=2;

b=[1];

for i=1:n-2

b=[b;0];

end

b=[b;1];

x0=(linspace(0,0,n))';

tol=10^(-6);

A;

b;

%下面采用算法8.2

for k=1:10000

x=[];

s=0;

for j=2:n

s=s+A(1,j)*x0(j);

end

x(1)=(b(1)-s)/A(1,1);

for i=2:n-1

s=0;t=0;

for j=1:i-1

s=s+A(i,j)*x(j);

end

for j=i+1:n

t=t+A(i,j)*x0(j);

end

x(i)=(b(i)-s-t)/A(i,i);

end

s=0;

for j=1:n-1

s=s+A(n,j);

end

x(n)=(b(n)-s)/A(n,n);

a=0;

switch(p)%根据判定标准算判定系数a case 1

for i=1:n

a=a+((x(i)-x0(i))^2);

end

a=sqrt(a);

case 2

xx=0;

for i=1:n

a=a+((x(i)-x0(i))^2);

xx=xx+x(i)^2;

end

a=sqrt(a)/sqrt(xx);

case 3

bb=A*(x)';

for i=1:n

a=a+((bb(i)-b(i))^2);

end

a=sqrt(a);

case 4

bb=A*(x)';

bbb=0;

for i=1:n

a=a+((bb(i)-b(i))^2);

bbb=bbb+bb(i)^2;

end

a=sqrt(a)/sqrt(bbb);

end

if a<=tol

x

k

break

else

x0=x;

end

end

SOR：

format long

n=input('请输入n的值：'); w=input('请输入松弛因子');

p=input('请选择判定标准代号：(1:相邻差量绝对误差；2：相邻差量相对误差；3：残量绝对误差；4：残量相对误差)');

A=zeros(n);

for i=1:n-1

A(i,i)=2;

A(i+1,i)=-1;

A(i,i+1)=-1;

end

A(n,n)=2;

b=[1];

for i=1:n-2

b=[b;0];

end

b=[b;1];

x0=(linspace(0,0,n))';

tol=10^(-6);

A;

b;

%下面采用算法8.3

for k=1:10000

x=[];

s=0;

for j=2:n

s=s+A(1,j)*x0(j);

end

x(1)=(1-w)*x0(1)+w*(b(1)-s)/A(1,1);

for i=2:n-1

s=0;t=0;

for j=1:i-1

s=s+A(i,j)*x(j);

end

for j=i+1:n

t=t+A(i,j)*x0(j);

end

x(i)=(1-w)*x0(i)+w*(b(i)-s-t)/A(i,i);

end

s=0;

for j=1:n-1