数学实验 5：线性代数方程组的数值解法

实验 5：线性代数方程组的数值解法

习题3：

已知方程组Ax b =,其中20*20

A R

∈，定义为：

3

1/21/41/231/21/41/41/231/21/41/41/231/21/4

1/23--⎡⎤⎢⎥---⎢

⎥⎢⎥---⎢⎥-⎢⎥⎢⎥

---⎢⎥

--⎣

⎦

试通过迭代法求解此方程组，认识迭代法收敛的含义以及迭代初值和方程组系数矩阵性质对

收敛速度的影响。实验要求：

（1） 选取不同的初始向量x0和不同的方程组右端向量b ，给定迭代误差要求，用雅可比

迭代法和高斯-赛德尔迭代法计算，观测得到的迭代向量序列是否均收敛？若收敛，记录迭代次数，分析计算结果并得出结论；

（2） 取定右端向量b 和初始向量x0，将A 的主对角线元素成倍的增长若干次，非主对角

元素不变，每次用雅可比迭代法计算，要求迭代误差满足(1)()5

||||10k k x x +-∞-<,比

较收敛速度，分析现象并得出结论。

1、 程序设计（可直接粘贴运行）

1） Jacobi 迭代法

function y=jacobi(a,b,x0,e,m)

%定义jacobi 函数，其中：a,b 为线性方程组Ax b =中的矩阵和右端向量；x0为初始值； %e 和m 分别为人为设定的精度和预计迭代次数；运行结果y 为迭代的结果和所有中间值组成的 %矩阵 y=0;

%对y 初始化

d=diag(diag(a)); %按雅可比迭代标准形形式取主对角元素作为矩阵D u=-triu(a,1); %取上三角矩阵u l=-tril(a,-1);

%取下三角矩阵l

bj=d^-1*(l+u); fj=d^-1*b;

x=[x0,zeros(20,m-1)];

%初始化x,其中x1=x0,即初始值

for k=1:m

%人为规定迭代次数，防止不收敛迭代导致死循环

x(:,k+1)=bj*x(:,k)+fj; %jacobi 迭代 if norm(x(:,k+1)-x(:,k),inf)

%判断迭代后是否满足迭代中止条件：(1)()||||k k x x e

+∞-<

y=x(:,1:k+1); %赋给y 所有中间值和迭代结果 sizej=k ; %若去掉;号，则输出迭代次数 break

%并结束迭代

end

%若不成立，继续迭代

end

%以下部分为验证迭代公式收敛的方法，仅需运行一次即可，因为收敛性完全由A 矩阵决定,而A %在本题是固定不变的；通过判断x Bx f =+中B 的谱半径或范数大小（B 在jacobi 迭代法中%

为矩阵bj ），即可得知收敛性:

e=eig(bj)

%输出全部特征值i λ，若max ||()1

i B λρ=<，则收敛

n1=norm(bj,1); %计算1-范数 n2=norm(bj);

%计算2-范数 nn=norm(bj,inf)；

%计算∞-范数 q=min([n1 n2 nn])

%由于谱半径不超过人以一种范数，所以只要范数的最小值q<1,也可判断迭代法收敛

2） Gauss 迭代法：与Jacobi 程序结构相同，不再注释

function y=gauss(a,b,x0,e,m) y=0;

d=diag(diag(a)); u=-triu(a,1); l=-tril(a,-1); x=[x0,zeros(20,m-1)]; bgs=(d-l)^-1*u; fgs=(d-l)^-1*b;

for k=1:m

x(:,k+1)=bgs*x(:,k)+fgs;

if norm(x(:,k+1)-x(:,k),inf)

y=x(:,1:k+1); sizeg=k; break end end

e=eig(bgs) n1=norm(bgs,1); n2=norm(bgs); nn=norm(bgs,inf); min([n1 n2 nn])

3） 操作函数：

%构造矩阵A

n=20;

a1=sparse(1:n,1:n,3,n,n); %按稀疏矩阵的输入法构造，比较方便 a2=sparse(1:n-1,2:n,-0.5,n,n); a3=sparse(1:n-2,3:n,-0.25,n,n); a=a1+a2+a3+a2'+a3'; a=full(a);

%还原为满矩阵

%通过给定不同的初始向量x0或者右端项b ，以及规定不同的误差要求，进行jacobi 和gauss %迭代，得到的结果y1、y2位两种迭代的次数，同时输出迭代结果，便于分析 b= x0= e= m=

y1=jacobi(a,b,x0,e,m); y2=gauss(a,b,x0,e,m);

4） 改变矩阵A ：

先对jacobi 函数作一定修正，方便分析，命名为jacobi2,如下:

function y=jacobi2(a,b,x0,e,m) d=diag(diag(a)); u=-triu(a,1); l=-tril(a,-1); bj=d^-1*(l+u); fj=d^-1*b;

x=[x0,zeros(20,m-1)];

n1=norm(bj,1);

%计算范数

n2=norm(bj); nn=norm(bj,inf); q=min([n1 n2 nn]);

y(1)=q;

%输出结果1：范数的最小值，判断收敛速度的方法

for k=1:m

x(:,k+1)=bj*x(:,k)+fj;

if norm(x(:,k+1)-x(:,k),inf)

%输出结果2：迭代次数

break end end

%改变A 矩阵主对角元素的值，比较jacobi 迭代的收敛速度，即迭代误差满足%(1)

()5||||10k k x

x +-∞-<时的迭代次数

b=(1:20)';

%取定b