Jacobi Gauss-Seidel 共轭梯度算法比较

(k )
=0 或 || r
(k )
|| ，则终止计算，输出“第 k 次迭代”以及 x* x( k ) ，否则转 3.
x=zeros(N,1);%迭代近似向量 eps=0.0000001;%精度 r=b-A*x;d=r; for k=0:N-1 fprintf('第%d次迭代',k+1); a=(norm(r)^2)/(d'*A*d) x=x+a*d rr=b-A*x; break; end C=(norm(rr)^2)/(norm(r)^2); d=rr+C*d; %rr=r(k+1) if (norm(rr)<=eps)||(k==N-1)
源程序：
clear;clc; A=[10,1,2,3,4; 1,9,-1,2,-3; 2,-1,7,4,-5; 3,2,3,12,-1; 4,-3,-5,-1,15]; b=[12;-27;14;-17;12]; N=length(b); fprintf('库函数计算结果: '); x=inv(A)*b x=[0,0,0,0,0]'; D=diag(diag(A)); E=-tril(A,-1); F=-triu(A,1); B=inv(D-E)*F; g=inv(D-E)*b; eps=0.001; for k=1:50 fprintf('第%d次迭代： ',k); y=B*x+g; fprintf('\n与上次计算结果距离（2 范数）：%f\n',norm(x-y)^2); if norm(x-y)<eps break; end x=y end
x=B*x0+f; i=i+1; end x i
运行结果：
（2） Gauss-Seidel 迭代法 Gauss-Seidel 迭代法算法：
(0) (0) (0) (0) 1. 输入 A (aij ), b (b1,…，bn ) ， x(0) ( x1 , x2 , x3 ,…，xn ) ，精度 ，维数 N.
题目二：用共轭梯度法或最速下降法求解 Hilbert 矩阵问题。 （程老师上课布置的作业）
解： 由于题目一种已使用共轭梯度法，在题目二中采用最速下降法求解。 算法： 按照每一个元素以及 Hilbert 矩阵的元素特性，即 amn
1 ，建立 Hilbert 矩阵，对其 m n 1
中元素用循环语句逐一建立。要想完成此题，我们需要先建立一个 Hilbert 矩阵，在此我以 5 ×5 的 Hilbert 矩阵为例，并令非齐次项矩阵 b=[1,1,1,1,1]’. 在快速下降算法中，设从 x
迭代初始向量取 x(0) (0,0,0,0,0)T .
解： （1） Jacobi 迭代法算法：
(0) (0) (0) (0) 1. 输入 A (aij ), b (b1,…，bn ) ， x(0) ( x1 , x2 , x3 ,…，xn ) ，精度 .
2. 置 i=1
B I D A， 3. 对 k=1,2,„， n， D diag (a11 , a22 ,…, ann ) ， g D1b ， x( k 1) Bx( k ) g
5. 若 r 6. 源程序：
clear A=[10,1,2,3,4; 1,9,-1,2,-3; 2,-1,7,4,-5; 3,2,3,12,-1; 4,-3,-5,-1,15]; b=[12;-27;14;-17;12];%自定义初始化非齐次部分 N=length(b); X=inv(A)*b %解向量的维数 %库函数计算结果 fprintf('库函数计算结果：');
function [x,n]=fastdown(A,b,x0,eps) %采用最速下降法求解关于线性Hilbert矩阵Ax=b的解 %Hilbert矩阵:A %线性方程组中非齐次项的向量:b %迭代初始向量:x0 %解的精度:eps %线性方程组的解:x %求出所需精度的解得实际迭代步数:n K=5; A=zeros(K,K); for m=1:K for n=1:K A(m,n)=1/(m+n-1); end end b=[1;1;1;1;1]; eps = 1.0e-6; x0=[0,0,0,0,0]'; r d x = b-A*x0; = dot(r,r)/dot(A*r,r); = x0+d*r;
(0)
出发，已迭代计算至 x
(k )
， ( x) 在 x
(k )
处的负梯度方向
，沿此方向 的函数值可表示为 r ( k ) ( x( k ) ) b Ax( k ) （ 称 为 搜 索 方 向 ）
( x( kΒιβλιοθήκη ) r k( ))，为求极小点令
d ( x ( k ) r ( k ) ) Ax( k ) b, r ( k ) Ar ( k ) b, r ( k ) 0 d
1
4. 若 || x x
(0)
|| ，输出最终结果 x 和迭代次数 i ，停机；否则转 5.
(0) 5. 置 i+1 i， xk xk （k=1, 2, …n） ，转 3；否则，停机
源程序：
clear clc i=1; a=[10,1,2,3,4; 1,9,-1,2,-3; 2,-1,7,4,-5; 3,2,3,12,-1; 4,-3,-5,-1,15]; d=diag(diag(a)); l=d-tril(a); u=d-triu(a); d0=inv(d); b=[12;-27;14;-17;12]; x0=[0,0,0,0,0]'; B=d0*(l+u); f=d0*b; x=B*x0+f; while norm(x-x0,inf)>=1e-4 x0=x;
r=rr; end
运算结果：
结果分析以及结论： 从结果来看，Jacobi 迭代法用了 45 次迭代得到精确解，Gauss-Seidel 迭代法用了 21 次迭代 得到精确解， 共轭梯度法仅仅用了 5 次就得到了精确解。 所以， 我们不难得出这样一个结论， 对于此题所需解决的线性方程组求解问题， 三种迭代方法都是收敛的， 且由迭代次数可以看 出， 共轭梯度法的收敛速度最快， 运用这种方法求解这道题的算法是最有效的， Gauss-Seidel 迭代法次之，Jacobi 迭代法最慢。
n=1; while(norm(x-x0)>eps) x0 = x; r d x n end = b-A*x0; = dot(r,r)/dot(A*r,r); = x0+d*r; = n + 1;
运行结果：
2. 置 K=1. 3. 由 Gauss-Seidel 算法， D diag (a11 , a22 ,…, ann )
0 a21 ������ =
0 0
0 0 0 … …
0 0 0 … an,n1
0 0 0 ,������ = … 0
0
0 0
a12
0
0
a13
a23 0 0
0
… … … …
得

(k )

r (k ) , r (k ) Ar ( k ) , r ( k )
于是 x( k 1) x( k ) k r ( k ) (k 1, 2,…) ，且序列 x 源程序：
满足 ( x
(k )
( k 1)
) ( x( k ) )(k 1, 2,…) .
第三章上机题目： 题目一：试分别用：
（ 1 ） Jacobi 迭 代 法 ； （ 2 ） Gauss-Seidel 迭 代 法 和 ； （3）共轭梯度法解线性方程组：
10 1 2 3 4
1 2 3 4 9 −1 2 −3 −1 7 3 −5 2 3 12 −1 −3 −5 −1 15
=
12 −27 14 −17 12
4. 对于 k=1,2,„,计算
r ( k ) b Ax( k )
k 1 r ( k ) , Ad ( k 1) / d ( k 1) , Ad ( k 1) , d ( k ) r ( k ) k 1d ( k 1)
k r ( k ) , d ( k ) / d ( k ) , Ad ( k ) , x ( k 1) x k k d ( k )
运行结果：
（3） 共轭梯度法 共轭梯度法算法：
(0) (0) (0) (0) 1. 输入 A (aij ), b (b1,…，bn ) ，初始向量 x(0) ( x1 , x2 , x3 ,…，xn ) ，误差限 .
2. 库函数计算结果作对比. 3. 计算 d
(0)
r (0) b Ax (0) , 0 r (0) , d (0) / d (0) , Ad (0) , x(1) x(0) 0 d (0) .
…
a1n
a2n , … an1,n
0
a31 … an1
a32 … an 2
0
0
0
0
x( k 1) ( D L)1Ux( k ) ( D L)1b
4. 若 || x x
(0)
|| ，输出“第 i 次迭代”和最终结果 x ，停机；否则转 5.
5. 置 k+1 k， xi xi(0) （i=1, 2, …n） ，转 3；否则，停机