通过为 qmr 提供用来计算 A*x 和 A'*x 的函数句柄(而非系数矩阵 A)来求解线性方程组。
创建一个非对称三对角矩阵。预览该矩阵。
A = gallery('wilk',21) + diag(ones(20,1),1)
A = 21×21
10 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
1 9 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 1 8 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 1 7 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 1 6 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 1 5 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 1 4 2 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 1 3 2 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 1 2 2 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 1 1 2 0 0 0 0 0 0 0 0 0 0
⋮
由于此三对角矩阵有特殊的结构,您可以用函数句柄来表示 A*x 运算。当 A 乘以向量时,所得向量中的大多数元素为零。结果中的非零元素对应于 A 的非零三对角元素。
表达式 Ax 变为:
Ax=[1020⋯⋯01920⋮01⋱20⋮010⋱⋱⋮0⋱1⋱0⋮⋱⋱⋱20⋯⋯0110][x1x2x3⋮⋮x21]=[10x1+2x2x1+9x2+2x3⋮⋮x19+9x20+2x21x20+10x21]。
结果向量可以写为三个向量的和:
Ax=[10x1+2x2x1+9x2+2x3⋮⋮x19+9x20+2x21x20+10x21]=[0x1x2⋮x20]+[10x19x2⋮9x2010x21]+2⋅[x2x3⋮x210]。
同样,ATx 的表达式变为:
ATx=[1010⋯⋯02910⋮02⋱10⋮020⋱⋱⋮0⋱1⋱0⋮⋱⋱⋱10⋯⋯0210][x1x2x3⋮⋮x21]=[10x1+x22x1+9x2+x3⋮⋮2x19+9x20+x212x20+10x21]。
ATx=[10x1+x22x1+9x2+x3⋮⋮2x19+9x20+x212x20+10x21]=2⋅[0x1x2⋮x20]+[10x19x2⋮9x2010x21]+[x2x3⋮x210]。
在 MATLAB® 中,编写一个函数来创建这些向量并将它们相加,从而根据标志输入给出 A*x 或 A'*x 的值:
function y = afun(x,flag)
if strcmp(flag,'notransp') % Compute A*x
y = [0; x(1:20)] ...
+ [(10:-1:0)'; (1:10)'].*x ...
+ 2*[x(2:end); 0];
elseif strcmp(flag,'transp') % Compute A'*x
y = 2*[0; x(1:20)] ...
+ [(10:-1:0)'; (1:10)'].*x ...
+ [x(2:end); 0];
end
end
(该函数作为局部函数保存在示例的末尾。)
现在,通过为 qmr 提供用于计算 A*x 和 A'*x 的函数句柄,求解线性方程组 Ax=b。使用容差 1e-6 和 25 次迭代。指定 b 为 A 的行总和,使得 x 的实际解是由 1 组成的向量。
b = full(sum(A,2));
tol = 1e-6;
maxit = 25;
x1 = qmr(@afun,b,tol,maxit)
qmr converged at iteration 19 to a solution with relative residual 4.7e-07.
x1 = 21×1
1.0000
1.0000
1.0000
1.0000
1.0000
1.0000
1.0000
1.0000
1.0000
1.0000
⋮
局部函数
function y = afun(x,flag)
if strcmp(flag,'notransp') % Compute A*x
y = [0; x(1:20)] ...
+ [(10:-1:0)'; (1:10)'].*x ...
+ 2*[x(2:end); 0];
elseif strcmp(flag,'transp') % Compute A'*x
y = 2*[0; x(1:20)] ...
+ [(10:-1:0)'; (1:10)'].*x ...
+ [x(2:end); 0];
end
end