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

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

习题3： 已知方程组Ax?b,其中A?R20*20，定义为：

?1/2?1/4?3??1/23?1/2???1/4?1/23????1/4???1/4?1/2?1/2?1/4??????1/4?3?1/2???1/23?

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

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

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

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

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

元素不变，每次用雅可比迭代法计算，要求迭代误差满足较收敛速度，分析现象并得出结论。

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

1） Jacobi迭代法

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

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

%对y初始化

%按雅可比迭代标准形形式取主对角元素作为矩阵D %取上三角矩阵u %取下三角矩阵l

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)];

for k=1:m

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

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

%判断迭代后是否满足迭代中止条件： y=x(:,1:k+1); sizej=k ; break

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

||x(k?1)?x(k)||??10?5,比

||x(k?1)?x(k)||??e

%赋给y所有中间值和迭代结果 %若去掉;号，则输出迭代次数 %并结束迭代

end end

%若不成立，继续迭代

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

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

e=eig(bj)

%输出全部特征值?i，若n1=norm(bj,1); n2=norm(bj);

max|?i|??(B)?1

，则收敛

%计算1-范数 %计算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;

for k=1:m

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

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

%输出结果2：迭代次数

%计算范数

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

%改变A矩阵主对角元素的值，比较jacobi迭代的收敛速度，即迭代误差满足%||x(k?1)?x(k)||??10?5时的迭代次数

%取定b

b=(1:20)';

x0=20*ones(20,1); e=10^-5; r=20; y=0; m=50; for k=1:r;

%取定x0 %取定精度

%设置主对角元素扩大的最大倍数

a1=k*sparse(1:n,1:n,3,n,n); %只需更改A的主对角元素 a=a1+a2+a3+a2'+a3';

a=full(a); end y

%输出对角线扩大倍数\\最小范数\\迭代次数

%重新构造A

y(k,1:3)=[k,jacobi2(a,b,x0,e,m)];

2、 运行结果分析

1） 不同初值（x0）和右端项（b）条件下的解的情况

表一：收敛性判断 Jacobi Gauss ?(B) 0.0163 0.0008 1-范数 0.0167 0.0084 2-范数 0.0163 0.0084 ?-范数 0.0167 0.0084 可以看到，矩阵A无论是谱半径或是任意范数的值都小于1，可知在A不变的情况下，Jacobi和Gauss法必然收敛。

2） b取不同的值，x0=20*ones(20,1), e=10^-5, m=50 条件下的情况对比 迭代次数 Jacobi Gauss B=[1:20]’ 24 16 B=[10:10:200]’ 26 17 B=[20:-1:1]’ 24 16 B=20*ones(20,1) 23 15 B=2000*ones(20,1) 30 20 根据1）分析的结果，可以证明无论b取任何值，采用两种方法迭代均收敛，但b的值的变化会影响迭代的次数；且Gauss迭代法总是比Jacobi迭代法收敛速度更快。 下表列出的是B=[1:20]’情况下部分结算结果，可以很明显的看到两种迭代法的收敛速度不同: Jacobi K=1 X1 X2 X3 Gauss X1 X2 X3 5.3333 9.0000 K=1 5.3333 6.5556 7.5370 K=2 K=3 K=4 K=5 K=6 … K=22 … 0.7247 … 1.3444 … 2.0072 … K=14 … 0.7247 … 1.3444 … 2.0072 K=23 0.7247 1.3444 2.0072 K=15 0.7247 1.3444 2.0072 K=24 0.7247 1.3444 2.0072 K=16 0.7247 1.3444 2.0072 标准值 0.7247 1.3444 2.0072 标准值 0.7247 1.3444 2.0072 2.7500 1.5394 1.0874 0.8858 0.7982 4.3333 2.6644 1.9248 1.6082 1.4652 K=2 K=3 K=4 K=5 K=6 11.0000 5.8056 3.7199 2.7801 2.3626 2.1717 2.0540 1.1354 0.8539 0.7657 0.7377 2.9432 1.8459 1.5033 1.3949 1.3605 3.7386 2.5553 2.1815 2.0627 2.0249

由于精度问题，在迭代的最后几次中从显示的数位已经不能看出标准值与计算值得差别，但是若采用long显示设定，就可以看到更多位小数的显示，其结果符合最初设定的精度e，数据繁琐，略。

3） x0取不同的值，b=20*ones(20,1), e=10^-5, m=50 条件下的情况对比 迭代次数 Jacobi Gauss X0=[1:20]’ 22 15 X0=[10:10:200]’ 27 18 X0=[20:-1:1]’ X0=20*ones(20,1) X0=2000*ones(20,1) 22 15 23 15 31 20 根据1）分析的结果，可以证明无论X0取任何值，采用两种方法迭代均收敛，但x0的值的变化会影响迭代的次数；且Gauss迭代法总是比Jacobi迭代法收敛速度更快。

下表列出的是x0=[1:20]’情况下部分结算结果，可以很明显的看到两种迭代法的收敛速 度不同: Jacobi X1 X2 X3 Gauss X1 X2 X3

4） b=(1:20)';x0=20*ones(20,1);e=10^-5; m=50;（固定各个参数）

r=20;（改变a的主对角元素，依次扩大1倍、2倍……20倍） 运行结果：

主对角元素扩大倍数 1.0000 2.0000 3.0000 4.0000 5.0000 6.0000 7.0000 8.0000 9.0000 10.0000 11.0000 12.0000 13.0000 14.0000 15.0000 16.0000 17.0000 18.0000 19.0000 20.0000

三个范数的最小值q 0.4893 0.2447 0.1631 0.1223 0.0979 0.0816 0.0699 0.0612 0.0544 0.0489 0.0445 0.0408 0.0376 0.0350 0.0326 0.0306 0.0288 0.0272 0.0258 0.0245 迭代次数 21.0000 12.0000 9.0000 8.0000 8.0000 7.0000 7.0000 7.0000 6.0000 6.0000 6.0000 6.0000 6.0000 6.0000 6.0000 6.0000 6.0000 6.0000 5.0000 5.0000 K=1 7.2500 7.6667 8.1667 K=1 7.2500 8.7083 9.8056 K=2 K=3 K=4 K=5 K=6 … … … … K=20 9.6327 K=21 9.6327 K=22 9.6327 标准值 9.6327 8.6250 9.2228 9.4536 9.5541 9.5974 … 9.9583 10.814 11.183 11.341 11.411 K=2 K=3 K=4 K=5 K=6 10.757 11.816 12.282 12.487 12.579 … 8.9352 9.4267 9.5720 9.6150 9.6276 … 10.654 11.228 11.398 11.448 11.463 11.813 12.408 12.584 12.636 12.650 … 11.4683 11.4683 11.4683 11.4683 12.6560 12.6560 12.6560 12.6560 K=13 9.6327 K=14 9.6327 K=15 9.6327 标准值 9.6327 11.4683 11.4683 11.4683 11.4683 12.6560 12.6560 12.6560 12.6560

本文来源：https://www.bwwdw.com/article/4hih.html