上海网站建设制作微信,数字营销软件,动易网站模板制作方法,外贸网站推广的方法目录
语法
说明
示例
线性系统的迭代解
使用指定了预条件子的 pcg
提供初始估计值
使用函数句柄代替数值矩阵 pcg函数的功能是求解线性系统 - 预条件共轭梯度法。
语法
x pcg(A,b)
x pcg(A,b,tol)
x pcg(A,b,tol,maxit)
x pcg(A,b,tol,maxit,M)
x pcg(A,b,tol,ma…目录
语法
说明
示例
线性系统的迭代解
使用指定了预条件子的 pcg
提供初始估计值
使用函数句柄代替数值矩阵 pcg函数的功能是求解线性系统 - 预条件共轭梯度法。
语法
x pcg(A,b)
x pcg(A,b,tol)
x pcg(A,b,tol,maxit)
x pcg(A,b,tol,maxit,M)
x pcg(A,b,tol,maxit,M1,M2)
x pcg(A,b,tol,maxit,M1,M2,x0)
[x,flag] pcg(___)
[x,flag,relres] pcg(___)
[x,flag,relres,iter] pcg(___)
[x,flag,relres,iter,resvec] pcg(___)
说明 x pcg(A,b) 尝试使用预条件共轭梯度法求解关于 x 的线性系统 A*x b。如果尝试成功pcg 会显示一条消息来确认收敛。如果 pcg 无法在达到最大迭代次数后收敛或出于任何原因暂停则会显示一条包含相对残差 norm(b-A*x)/norm(b) 以及该方法停止时的迭代次数的诊断消息。 x pcg(A,b,tol) 指定该方法的容差。默认容差是 1e-6。 x pcg(A,b,tol,maxit) 指定要使用的最大迭代次数。如果 pcg 无法在 maxit 次迭代内收敛将显示诊断消息。 x pcg(A,b,tol,maxit,M) 指定预条件子矩阵 M 并通过有效求解关于 y 的方程组 来计算 x其中 且 。该算法不显式形成 H。使用预条件子矩阵可以改善问题的数值属性和计算的效率。 x pcg(A,b,tol,maxit,M1,M2) 指定预条件子矩阵 M 的因子使得 M M1*M2。 x pcg(A,b,tol,maxit,M1,M2,x0) 指定解向量 x 的初始估计值。默认值为由零组成的向量。 [x,flag] pcg(___) 返回一个标志指示算法是否成功收敛。当 flag 0 时收敛成功。您可以将此输出语法用于之前的任何输入参数组合。如果指定了 flag 输出pcg 将不会显示任何诊断消息。 [x,flag,relres] pcg(___) 还会返回相对残差 norm(b-A*x)/norm(b)。如果 flag 为 0则 relres tol。 [x,flag,relres,iter] pcg(___) 还会返回计算出 x 时的迭代次数 iter。 [x,flag,relres,iter,resvec] pcg(___) 还会在每次迭代中返回残差范数向量包括第一个残差 norm(b-A*x0)。
示例
线性系统的迭代解 使用采用默认设置的 pcg 求解系数矩阵为方阵的线性系统然后在求解过程中调整使用的容差和迭代次数。 创建一个随机对称稀疏矩阵 A。还要为 Axb 的右侧创建一个由 A 的行总和组成的向量 b使得实际解 x 是由 1 组成的向量。
rng default
A sprand(400,400,.5);
A A*A;
b sum(A,2); 使用 pcg 求解 Axb。输出显示包括相对残差 ‖b−Ax‖/‖b‖ 的值。
x pcg(A,b);
pcg stopped at iteration 20 without converging to the desired tolerance 1e-06
because the maximum number of iterations was reached.
The iterate returned (number 20) has relative residual 3.6e-06. 默认情况下pcg 使用 20 次迭代和容差 1e-6对于此矩阵算法无法在 20 次迭代后收敛。然而残差接近容差因此算法可能只需更多迭代即可收敛。 使用容差 1e-7 和 150 次迭代再次求解方程组。
x pcg(A,b,1e-7,150);
pcg converged at iteration 129 to a solution with relative residual 1e-07.
使用指定了预条件子的 pcg 检查使用指定了预条件子矩阵的 pcg 来求解线性系统的效果。 创建一个对称正定带状系数矩阵。
A delsq(numgrid(S,102)); 定义 b 作为线性方程 Axb 的右侧。
b ones(size(A,1),1);
设置容差和最大迭代次数。
tol 1e-8;
maxit 100;
使用 pcg 根据请求的容差和迭代次数求解。指定五个输出以返回有关求解过程的信息 x 是计算 A*x b 所得的解。 fl0 是指示算法是否收敛的标志。 rr0 是计算的解 x 的相对残差。 it0 是计算 x 时所用的迭代次数。 rv0 是 ‖b−Ax‖ 的残差历史记录组成的向量。
[x,fl0,rr0,it0,rv0] pcg(A,b,tol,maxit);
fl0
fl0 1
rr0
rr0 0.0131
it0
it0 100 由于 pcg 未在请求的 100 次迭代内收敛至请求的容差 1e-8因此 fl0 为 1。 为了有助于缓慢收敛可以指定预条件子矩阵。由于 A 是对称矩阵请使用 ichol 生成预条件子 。通过指定 L 和 L 作为 pcg 的输入求解预条件方程组。
L ichol(A);
[x1,fl1,rr1,it1,rv1] pcg(A,b,tol,maxit,L,L);
fl1
fl1 0
rr1
rr1 8.0992e-09
it1
it1 79 在第 79 次迭代中使用 ichol 预条件子产生的相对残差小于规定的容差 1e-8 。输出 rv1(1) 是 norm(b) 且 rv1(end) 是 norm(b-A*x1)。 现在使用 michol 选项创建修正的不完全 Cholesky 预条件子。
L ichol(A,struct(michol,on));
[x2,fl2,rr2,it2,rv2] pcg(A,b,tol,maxit,L,L);
fl2
fl2 0
rr2
rr2 9.9618e-09
it2
it2 47 该预条件子优于使用零填充不完全 Cholesky 分解为本例中的系数矩阵生成的预条件子因此 pcg 能够更快地收敛。 通过绘制从初始估计值迭代编号 0开始的每次残差历史记录可以查看预条件子如何影响 pcg 的收敛速度。为指定的容差添加一个线条。
semilogy(0:length(rv0)-1,rv0/norm(b),-o)
hold on
semilogy(0:length(rv1)-1,rv1/norm(b),-o)
semilogy(0:length(rv2)-1,rv2/norm(b),-o)
yline(tol,r--);
legend(No Preconditioner,Default ICHOL,Modified ICHOL,Tolerance,Location,East)
xlabel(Iteration number)
ylabel(Relative residual)
如图所示 提供初始估计值 检查向 pcg 提供解的初始估计值的效果。 创建一个三对角稀疏矩阵。使用每行的总和作为 Axb 右侧的向量使 x 的预期解是由 1 组成的向量。
n 900;
e ones(n,1);
A spdiags([e 2*e e],-1:1,n,n);
b sum(A,2); 使用 pcg 求解 Axb 两次一次是使用默认的初始估计值一次是使用解的良好初始估计值。对两次求解均使用 200 次迭代和默认容差。将第二种求解中的初始估计值指定为所有元素都等于 0.99 的向量。
maxit 200;
x1 pcg(A,b,[],maxit);
pcg converged at iteration 35 to a solution with relative residual 9.5e-07.
x0 0.99*e;
x2 pcg(A,b,[],maxit,[],[],x0);
pcg converged at iteration 7 to a solution with relative residual 8.7e-07. 在这种情况下提供初始估计值可以使 pcg 更快地收敛。
返回中间结果 还可以通过在 for 循环中调用 pcg 来使用初始估计值获得中间结果。每次调用求解器都会执行几次迭代并存储计算出的解。然后将该解用作下一批迭代的初始向量。 例如以下代码会循环执行四次每次执行 100 次迭代并在 for 循环中每通过一次后均存储解向量
x0 zeros(size(A,2),1);
tol 1e-8;
maxit 100;
for k 1:4[x,flag,relres] pcg(A,b,tol,maxit,[],[],x0);X(:,k) x;R(k) relres;x0 x;
end X(:,k) 是在 for 循环的第 k 次迭代时计算的解向量R(k) 是该解的相对残差。
使用函数句柄代替数值矩阵 通过为 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×2020 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 01 19 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 00 1 18 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 00 0 1 17 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 00 0 0 1 16 1 0 0 0 0 0 0 0 0 0 0 0 0 0 00 0 0 0 1 15 1 0 0 0 0 0 0 0 0 0 0 0 0 00 0 0 0 0 1 14 1 0 0 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 1 13 1 0 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 1 12 1 0 0 0 0 0 0 0 0 0 00 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 变为 结果向量可以写为三个向量的和 在 MATLAB® 中编写一个函数来创建这些向量并将它们相加从而给出 A*x 的值
function y afun(x)
y [0; x(1:19)] ...[(20:-1:1)].*x ...[x(2:20); 0];
end 该函数作为局部函数保存在示例的末尾。 现在通过为 pcg 提供用于计算 A*x 的函数句柄求解线性系统 Axb。使用容差 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×10.04760.04750.05000.05260.05550.05880.06250.06660.07140.0769⋮ 检查 afun(x1) 是否产生由 1 组成的向量。
afun(x1)
ans 20×11.00001.00001.00001.00001.00001.00001.00001.00001.00001.0000⋮
局部函数
function y afun(x)
y [0; x(1:19)] ...[(20:-1:1)].*x ...[x(2:20); 0];
end
预条件共轭梯度法 预条件共轭梯度法 (PCG) 旨在利用对称正定矩阵的结构。其他一些算法也可以对对称正定矩阵执行运算但 PCG 在求解这些类型的方程组方面是最快和最可靠的 [1]。
提示 大多数迭代方法的收敛取决于系数矩阵的条件数 cond(A)。当 A 是方阵时可以使用 equilibrate 来改进其条件数它本身就能使大多数迭代求解器更容易收敛。但如果您随后会对经平衡处理的矩阵 B R*P*A*C 进行因式分解使用 equilibrate 还可以获得质量更好的预条件子矩阵。 可以使用矩阵重新排序函数如 dissect 和 symrcm来置换系数矩阵的行和列并在系数矩阵被分解以生成预条件子时最小化非零值的数量。这可以减少后续求解预条件线性系统所需的内存和时间。
参考 [1] Barrett, R., M. Berry, T. F. Chan, et al., Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods, SIAM, Philadelphia, 1994.
