数值代数上机习题：线性方程组的古典迭代解法

x ( k ) BGS x ( k 1) g GS ,
其分量形式为：
k 1,2,3,...
(k ) 1 n ( k 1) b1 ), ( a1 j x j x1 a j 2 11 n ( k ) 1 i 1 (k ) ( k 1) bi ), x ( a x aij x j i ij j a j j i 1 1 ii (k ) 1 n 1 ( k 1) bn ), xn ( anj x j a j 1 nn
4.SOR 迭代中最佳松弛因子的选取 由于 SOR 迭代法的效果和其松弛因子 w 的选取有关， 所以有必要选取合适的松弛因子。 《数 值线性代数》 （徐树方，北京大学出版社）中的相关结论，当选择最佳松弛因子
wopt
2 1 1 ( B) 2
时，SOR 方法的迭代速度最快。
算法使用 matlab 实现， 程序中采用无穷范数 || . || ， 每种迭代法最大允许的迭代次数为 100000 次。
b(1) ah 2 y (0) * b(n 1) ah 2 y (1) * ( h)
3.迭代终止条件 首先确定要求的精度 tol ，我们希望当
|| xk x* || tol
则停止迭代。 由《数值线性代数》 （徐树方，北京大学出版社） P105 定理 4.2.3，我们知道，对于迭代格 式 xk Bxk 1 g ，若 || B || q 1 且 || I || 1 ，则迭代序列 {xk } 的第 k 次近似解和精确解之 间有估计式 || xk x* ||
1 , n

yi 1 2 yi yi 1 yi 1 yi a h2 h
简化为：
( h) yi 1 (2 h) yi yi 1 ah 2
从而离散后得到的线性方程组的系数矩阵为：
h (2 h) (2 h) h A (2 h) h (2 h)
% ---------------------------- % % filename : jacobi.m % abstract : Jacobi迭代法解线性方程组 Ax = b % 调用方式: % [x, k, e] = jacobi(A,b,y) % 返回值： % x : 迭代解 % k : 迭代次数 % e : 与精确解的误差 % ---------------------------- % % Reference : 数值线性代数，徐树方，北京大学出版社 % ---------------------------- % function [x, k, e] = jacobi(A,b,y) N = 100000; tol = 1e-4; q D L U B g = = = = = = 1 - 1/(100*100); diag(diag(A)); triu(A)-A; tril(A)-A; D\(L+U); D\b; % y为精确解
然后在矩阵乘法时对下标处理一下即可。但是考虑到三种迭代方法的一般性，且本题中
n 100 并不是很大，所以实验中并没有采用紧缩存储，而是采用了直接存储。
2.边值条件的处理 由于差分得到的方程组的第一行和最后一行中分别出现了边值 y(0)与 y(1)作为常数项， 因此 要在常向量的第一项和最后一项作一些修改：
283 199 111 108
五．结果分析
1.迭代法的收敛速度比较
从表格数据可以看出，当 不是很小时，GS 迭代的速度要快于 Jacobi 迭代的速度，而 SOR 迭代（选取最佳松弛因子）的速度要远远快于前两者。当 减小时，三种迭代方法的收敛速 度均增快，而当 0.0001 时，求得 SOR 迭代的最佳松弛因子接近 1，基本退化为 GS 迭 代，这是因为当 很小时，方程系数矩阵的谱半径 ( A) 也很小，所以求出的最佳松弛因子
四．计算结果
题目要求 a
1 , n 100, 下面是 1,0.1,0.01,0.0001 时三种方法计算所得近似解需要的迭 2
代次数以及与精确解的误差。 1.Jacobi 迭代

迭代次数 k
与精确解的误差 || xk x* || 3.100e-4 8.826e-3 6.606e-2 4.950e-3
它的分量形式为：
xi
(k )
(1 ) xi
( k 1)

n 1 i 1 (k ) (k ) ( aij x j aij x j bi ), aii j 1 j i 1
其中 称为松弛因子，当 1 时称为超松弛；当 1 时叫低松弛； 1 时就是 G-S 迭 代。
wopt
2 1 1 ( B ) 2 接近于 1。
2.误差比较 从表格可以看出， 三种迭代方法所得近似解的误差基本一样。 当 1 时， 误差在 1e-4 左右， 符合定理 4.2.3 的结果。但当 减小时，误差却出现先增大再减小的情况，这可能是由于前
面确定迭代终止时的近似 q 1 h 不够准确所致。
3.超松弛迭代法（SOR 迭代法） ：
i 2,..., n 1
SOR迭代法是对G-S迭代法的加权平均改造，即 x
(k )
(1 ) x ( k 1) x ,
(k )
x
(k )
为G-S迭代解，即
x ( k ) (1 ) x ( k 1) ( D 1 Lx ( k ) D 1Ux ( k 1) D 1b),
2
参考文献 [1] 徐树方， 《数值线性代数》 ，北京大学出版社
附录：matlab 代码
% ---------------------------- % % filename : test.m % abstract : 用于测试的m脚本文件，输入系数eps， 显示三种迭代方法得到的近似解，迭代次数及与精确解的误差 % % ---------------------------- % % ---------------------------- % a = 0.5; n = 100; eps = input('请输入系数eps:'); [A b] = genmatrix(eps, a, n); y = exactsolution(eps, a, n); % Jacobi迭代 [x k e] = jacobi(A,b,y); disp('% Jacobi迭代 '); disp('近似解:'); disp(x); fprintf('误差：%.6f;迭代次数：%d\n\n',e,k); % GS迭代 % 输入eps % 生成系数矩阵和常向量 % 获得精确解
a13 a1n a23 a2 n 0 , an 1,n 0
则原方程组可改写为：
x BJ x g J ，
其中
（2.2）
BJ D 1 ( L U ), g J D 1b,
给定初始向量：
x0 ( x1 , x2 ,..., xn ) T ,
数值代数上机习题一：线性方程组的古典迭代解法
一．题目要求
考虑两点边值问题
d 2 y dy a ,0 a 1 2 dx dx y (0) 0, y (1) 1.
容易知道它的精确解为
y
1 a (1 e ) ax 1 / 1 e
x
为了把微分方程离散，把[0, 1]区间 n 等分，令 h 得到差分方程为：
q || xk 1 xk || 。 1 q
4
由 题 目 要 求 知 我 们 需 要 有 || xk x* || 10
，而由上面的迭代估计，只要
q 1 q 即可。而本题中 q 可近似取为 || xk 1 xk || 10 4 ，即 || xk 1 xk || 10 4 1 q q q 1 h 2 ，因此最后令迭代终止条件为 || xk 1 xk || 10 4 1 q q

1 0 .1 0.01 0.0001
22274 8754 662 122
2.GS 迭代

迭代次数 k
与精确解的误差 || xk x* || 2.900e-4 8.821e-3 6.606e-2 4.950e-3
1 0 .1 0.01 0.0001
( 0)
(0)
(0)
三．算法实现
1.方程的表示及存储
由于本题中线性方程组的系数矩阵为三对角矩阵，所以可以采用紧缩方法存储，即
0 A
(2 h) h (2 h) h (2 h) h (2 h) 0
D diag (a11 , a22 ,..., ann ),
0 a 21 L a31 an1 , 0
（2.1）
0 a32 0
an 2 an ,n 1
0 a12 0 U
要求： （1）对 1, a
1 , n 100, 分别用 Jacobi, G-S 和 SOR 迭代方法求线性方程组的解，要求 2
有 4 位有效数字，然后比较与精确解的误差。 （2）对 0.1, 0.01, 0.0001, 考虑同样的问题。
二．算法原理
1.Jacobi 迭代法（J-迭代法） ： 对于非奇异线性方程组 Ax b ，令 A D 1
2.Gauss-Seidel 迭代法（G-S 迭代法） ： 类似于 J-迭代法，给定初值：
x0 ( x1 , x2 ,..., xn ) T ,
( 0)
(0)
(0)
令
BGS ( D L) 1U , g GS ( D L) 1 b,
则得到 G-S 迭代公式：
从上面的介绍可以看出， 三种迭代法的本质区别在于迭代公式的不同， 因此算法的基本框架 如下：
1）输入初值 x0 ( x1 , x2 ,..., xn ) T , 精度要求 tol，最大允许迭代次数 N，令 k := 0,转 2); 2）由迭代公式计算出 x ( k 1) , 转 3); 3）若 || x ( k ) x ( k 1) || tol ,则停止，输出 x ( k 1) 作为方程的解；若否，转 4) 4）若 k = N，则停止，输出迭代失败信息；若否，令 k := k+1，转 2)
[x k e] = gs(A,b,y); disp('% GS迭代 '); disp('近似解:'); disp(x); fprintf('误差：%.6f;迭代次数：%d\n\n',e,k); % SOR迭代 [x k e] = sor(A,b,y); disp('% SOR迭代 '); disp('近似解:'); disp(x); fprintf('误差：%.6f;迭代次数：%d\n\n',e,k);
10788 3855 344 109
3.SOR 迭代

迭代次数 k
误差 || xk x* || 3.000e-4 8.829e-3 6.606e-2 4.950e-3
松 弛 因 子w 1.9384 1.8921 1.4985 1.0021
1 0 .1 0.01 0.0001
由（2.2）可以构造迭代公式：
( 0)
(0)
(0)
xk BJ xk 1 g ,
其分量形式为：
k 1,2,3,...
(k ) 1 n ( k 1) b1 ), ( a1 j x j x1 a11 j 2 n ( k ) 1 i 1 ( k 1) ( k 1) aij x j bi ), xi ( aij x j aii j 1 j i 1 (k ) 1 n 1 ( k 1) bn ), xn ( anj x j ann j 1