48．ichol函数

在MATLAB中，提供了ichol函数实现不完全Cholesky分解。函数的语法格式为：

L＝ichol（A）：使用零填充对A执行不完全Cholesky分解。

L＝ichol（A，opts）：使用opts指定的选项对A执行不完全Cholesky分解。

默认情况下，ichol引用A（稀疏矩阵）的下三角并生成下三角因子，其取值见表1-3。

表1-3　opts的取值及含义

【例1-49】本示例演示如何结合pcg使用预设子系统矩阵。

     %创建一个输入矩阵并尝试使用pcg求解方程组
     A = delsq(numgrid('S',100));
     b = ones(size(A,1),1);
     [x0,fl0,rr0,it0,rv0] = pcg(A,b,1e-8,100);

由于pcg未在请求的最多100次迭代内收敛至请求的容差1e－8，因此fl0为1。预设子条件可以使方程组更快地收敛。

     %使用只带一个输入参数的ichol构造一个含有零填充值的不完全Cholesky分解
     L = ichol(A);
     [x1,fl1,rr1,it1,rv1] = pcg(A,b,1e-8,100,L,L');

当使用零填充不完全Cholesky分解进行预调节时，由于pcg使相对残差在第77次迭代（it1的值）时趋于9.8e－09（rr1的值），且该值小于请求的容差1e－8，因此fl1为0。rv1（1）＝norm（b），rv1（78）＝norm（b－A∗x1）。

前一矩阵表示基于100×100网格、带Dirichlet边界条件的拉普拉斯算子离散化。这意味着使用修正的不完全Cholesky预设子条件可能获得更好效果。

     %使用michol选项创建修正的不完全Cholesky预设子条件
     L = ichol(A,struct('michol','on'));
     [x2,fl2,rr2,it2,rv2] = pcg(A,b,1e-8,100,L,L');

在这种情况下，只需要进行47次迭代就能获得收敛。

     %通过绘制从初始估计值(迭代数0)开始的每次残差历史记录
     %可以查看预设子条件如何影响pcg的收敛速度
     figure;
     semilogy(0:it0,rv0/norm(b),'b.');
     hold on;
     semilogy(0:it1,rv1/norm(b),'r.');
     semilogy(0:it2,rv2/norm(b),'k.');
     legend('没有预设子条件','IC(0)','MIC(0)');
     xlabel('迭代次数');
     ylabel('相对残差');
     hold off;

运行程序，效果如图1-20所示。

图1-20　预设子条件影响pcg收敛速度效果

彩色图片