matlab ichol,求解线性方程组 - 预条件共轭梯度法

通过为 pcg 提供用来计算 A*x 的函数句柄(而非系数矩阵 A)来求解线性方程组。

使用 gallery 生成一个 20×20 正定三对角矩阵。上对角线和下对角线上的元素都是 1，而主对角线上的元素是从 20 递减到 1。预览该矩阵。

n = 20;

A = gallery('tridiag',ones(n-1,1),n:-1:1,ones(n-1,1));

full(A)

ans = 20×20

20 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

1 19 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 1 18 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 1 17 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 1 16 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 1 15 1 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 1 14 1 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 1 13 1 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 1 12 1 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 1 11 1 0 0 0 0 0 0 0 0 0

⋮

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

表达式 Ax 变为：

[2010⋯⋯⋯00119100011810⋮⋮011710011610⋮⋮011510011410⋮⋮0113⋱000⋱⋱100⋯⋯⋯011][x1x2x3x4x5⋮⋮x20]=[20x1+x2x1+19x2+x3x2+18x3+x4⋮x18+2x19+x20x19+x20]。

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

[20x1+x2x1+19x2+x3x2+18x3+x4⋮x18+2x19+x20x19+x20]=[0x1⋮x19]+[20x119x2⋮x20]+[x2⋮x200]。

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

function y = afun(x)

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

[(20:-1:1)'].*x + ...

[x(2:20); 0];

end

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

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

b = ones(20,1);

tol = 1e-12;

maxit = 50;

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

pcg converged at iteration 20 to a solution with relative residual 4.4e-16.

x1 = 20×1

0.0476

0.0475

0.0500

0.0526

0.0555

0.0588

0.0625

0.0666

0.0714

0.0769

⋮

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

afun(x1)

ans = 20×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:19)] + ...

[(20:-1:1)'].*x + ...

[x(2:20); 0];

end