matlab进行稳定性分析,求解线性方程组 - 稳定双共轭梯度法

通过为 bicgstab 提供函数句柄来求解线性方程组，用函数句柄代替系数矩阵 A 来计算 A*x。

gallery 生成的 Wilkinson 测试矩阵之一是 21×21 三对角矩阵。预览该矩阵。

A = gallery('wilk',21)

A = 21×21

10 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

1 9 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 1 8 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 1 7 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 1 6 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 1 5 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 1 4 1 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 1 3 1 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 1 2 1 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 1 1 1 0 0 0 0 0 0 0 0 0 0

⋮

Wilkinson 矩阵有特殊的结构，因此您可以用函数句柄来表示 A*x 运算。当 A 乘以向量时，所得向量中的大多数元素为零。结果中的非零元素对应于 A 的非零三对角元素。此外，只有主对角线具有不等于 1 的非零值。

表达式 Ax 变为：

Ax=[1010⋯⋯⋯001910001810⋮⋮0171001610⋮⋮0151001410⋮⋮013⋱000⋱⋱100⋯⋯⋯0110][x1x2x3x4x5⋮⋮x21]=[10x1+x2x1+9x2+x3x2+8x3+x4⋮x19+9x20+x21x20+10x21]。

结果向量可以写为三个向量的和：

Ax=[0+10x1+x2x1+9x2+x3x2+8x3+x4⋮x19+9x20+x21x20+10x21+0]=[0x1⋮x20]+[10x19x2⋮10x21]+[x2⋮x210]。

在 MATLAB® 中，编写一个函数来创建这些向量并将它们相加，从而给出 A*x 的值：

function y = afun(x)

y = [0; x(1:20)] + ...

[(10:-1:0)'; (1:10)'].*x + ...

[x(2:21); 0];

end

(该函数作为局部函数保存在示例的末尾。)

现在，通过为 bicgstab 提供用于计算 A*x 的函数句柄，求解线性方程组 Ax=b。使用容差 1e-12 和 50 次迭代。

b = ones(21,1);

tol = 1e-12;

maxit = 50;

x1 = bicgstab(@afun,b,tol,maxit)

bicgstab converged at iteration 11.5 to a solution with relative residual 5.2e-13.

x1 = 21×1

0.0910

0.0899

0.0999

0.1109

0.1241

0.1443

0.1544

0.2383

0.1309

0.5000

⋮

检查 afun(x1) 是否产生由 1 组成的向量。

afun(x1)

ans = 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)

y = [0; x(1:20)] + ...

[(10:-1:0)'; (1:10)'].*x + ...

[x(2:21); 0];

end