vue.js和vs做网站比较,航拍类wordpress模板,idea制作网站,做网站时新闻的背景图目录
语法
说明
示例
不同类型的不完全 LU 分解
不完全 LU 分解的调降容差
使用 ilu 作为预条件子来求解线性系统 ilu函数的功能是对矩阵进行不完全 LU 分解。
语法
[L,U] ilu(A)
[L,U,P] ilu(A)
W ilu(A)
[___] ilu(A,options)
说明
[L,U] ilu(A) 用零填充执行…目录
语法
说明
示例
不同类型的不完全 LU 分解
不完全 LU 分解的调降容差
使用 ilu 作为预条件子来求解线性系统 ilu函数的功能是对矩阵进行不完全 LU 分解。
语法
[L,U] ilu(A)
[L,U,P] ilu(A)
W ilu(A)
[___] ilu(A,options)
说明
[L,U] ilu(A) 用零填充执行稀疏矩阵 A 的不完全 LU 分解并返回下三角矩阵 L 和上三角矩阵 U。
[L,U,P] ilu(A) 还返回置换矩阵 P并满足 L 和 U 是 P*A 或 A*P 的不完全因子。默认情况下P 是用于不使用主元消去的不完全 LU 分解的单位矩阵。
W ilu(A) 返回 LU 因子的非零值。输出 W 等于 L U - speye(size(A))。
[___] ilu(A,options) 使用结构体 options 指定的选项对 A 执行不完全 LU 分解。 例如通过将 options 的 type 字段设置为 ilutp您可以使用主元消去执行不完全 LU 分解。然后通过将 milu 字段设置为 row 或 col可以指定保留修正不完全 LU 分解的行值总和或列值总和。这种选项组合返回置换矩阵 P使得 L 和 U 是 row 选项的 A*P 的不完全因子其中 U 是列置换的或 L 和 U 是 col 的 P*A 的不完全因子其中 L 是行置换的。
[L,exitflag] logm(A) 返回描述 logm 的退出条件的标量 exitflag
示例
不同类型的不完全 LU 分解 ilu 函数提供三种类型的不完全 LU 分解零填充分解 (ILU(0))、Crout 版本 (ILUC)以及具有阈值调降和主元消去的分解 (ILUTP)。 默认情况下ilu 执行稀疏矩阵输入的零填充不完全 LU 分解。例如求具有 7840 个非零值的稀疏矩阵的完全和不完全分解。其完全 LU 因子有 126,478 个非零值其具有零填充的不完全 LU 因子有 7840 个非零值与 A 的数量相同。
A gallery(neumann,1600) speye(1600);
n nnz(A)n 7840n nnz(lu(A))n 126478n nnz(ilu(A))
n 7840 由于零填充分解能在其 LU 因子中保留输入矩阵的稀疏模式因此分解的相对误差在 A 的非零元素模式中实质上为零。
[L,U] ilu(A);
e norm(A-(L*U).*spones(A),fro)/norm(A,fro)
e 4.8874e-17 然而这些零填充因子的乘积并非原始矩阵的良好逼近。
e norm(A-L*U,fro)/norm(A,fro)
e 0.0601 为了提高精确度可以使用其他类型的具有填充的不完全 LU 分解。例如使用具有 1e-6 调降容差的 Crout 版本。
options.droptol 1e-6;
options.type crout;
[L,U] ilu(A,options); 不完全分解的 Crout 版本在其 LU 因子中具有 51,482 个非零值少于完全分解。在具有填充的情况下不完全 LU 因子的乘积是原始矩阵的更好逼近。
n nnz(ilu(A,options))
n 51482
e norm(A-L*U,fro)./norm(A,fro)
e 9.3040e-07 作为比较对于相同的输入矩阵 A具有阈值调降和主元消去的不完全分解将给出类似于 Crout 版本的结果。
options.type ilutp;
[L,U,P] ilu(A,options);
n nnz(ilu(A,options))
n 51541
norm(P*A-L*U,fro)./norm(A,fro)
ans 9.4960e-07
不完全 LU 分解的调降容差 更改不完全 LU 分解的调降容差以分解稀疏矩阵。 加载 west0479 矩阵它是一个非对称的 479×479 实数值稀疏矩阵。使用 condest 估计该矩阵的条件数。
load west0479
A west0479;
c1 condest(A)
c1 1.4244e12 使用 equilibrate 改进矩阵的条件数。对原始矩阵 A 进行置换和重新缩放以创建一个新矩阵 B R*P*A*C它具有更好的条件数且对角线元只有 1 和 -1。
[P,R,C] equilibrate(A);
B R*P*A*C;
c2 condest(B)
c2 5.1036e04 指定选项以执行带有阈值调降和主元消去的 B 的不完全 LU 分解保留行值总和不变。为了便于比较首先将填充的调降容差指定为零这将产生完全 LU 分解。
options.type ilutp;
options.milu row;
options.droptol 0;
[L,U,P] ilu(B,options); 这种分解在逼近输入矩阵 B 时非常准确但因子明显比 B 更稠密。
e norm(B*P-L*U,fro)/norm(B,fro)
e 1.0345e-16
nLU nnz(L)nnz(U)-size(B,1)
nLU 19118
nB nnz(B)
nB 1887 可以通过更改调降容差以获得不完全 LU 因子这些因子更稀疏但在逼近 B 时不太准确。例如以下绘图显示了不完全因子的密度与输入矩阵的密度之比以及不完全分解的相对误差它们分别相对于调降容差的图。
ntols 20;
tau logspace(-6,-2,ntols);
e zeros(1,ntols);
nLU zeros(1,ntols);
for k 1:ntolsoptions.droptol tau(k);[L,U,P] ilu(B,options);nLU(k) nnz(L)nnz(U)-size(B,1);e(k) norm(B*P-L*U,fro)/norm(B,fro);
end
figure
semilogx(tau,nLU./nB,LineWidth2)
title(Ratio of nonzeros of LU factors with respect to B)
xlabel(drop tolerance)
ylabel(nnz(L)nnz(U)-size(B,1)/nnz(B),FontNameFixedWidth)
如图所示 figure
loglog(tau,e,LineWidth2)
title(Relative error of the incomplete LU factorization)
xlabel(drop tolerance)
ylabel($||BP-LU||_F\,/\,||B||_F$,Interpreterlatex)
如图所示 在此示例中具有阈值调降的不完全 LU 分解的相对误差与调降容差处于相同的数量级但不能保证一定会发生此情况。
使用 ilu 作为预条件子来求解线性系统 使用不完全 LU 分解作为 bicgstab 的预条件子来求解线性系统。 加载 west0479它是一个非对称的 479×479 实稀疏矩阵。
load west0479
A west0479; 定义b以使 Axb 的实际解是全为 1 的向量。
b full(sum(A,2));
设置容差和最大迭代次数。
tol 1e-12;
maxit 20;
使用 bicgstab 根据请求的容差和迭代次数求解。指定五个输出以返回有关求解过程的信息 x 是计算 A*x b 所得的解。 fl0 是指示算法是否收敛的标志。 rr0 是计算的解 x 的相对残差。 it0 是计算 x 时所用的迭代次数。 rv0 是 ‖b−Ax‖ 的由每个二分之一迭代的残差历史记录组成的向量。
[x,fl0,rr0,it0,rv0] bicgstab(A,b,tol,maxit);
fl0
fl0 1
rr0
rr0 1
it0
it0 0 bicgstab 未在请求的 20 次迭代内收敛至请求的容差 1e-12因此 fl0 为 1。实际上bicgstab 的行为太差因此初始估计值 x0 zeros(size(A,2),1) 是最佳解并会返回最佳解如 it0 0 所示。 为了有助于缓慢收敛可以指定预条件子矩阵。由于 A 是非对称的请使用 ilu 生成预条件子 ML U。指定调降容差以忽略值小于 1e-6 的非对角线元。通过指定 L 和 U 作为 bicgstab 的输入求解预条件方程组
options struct(type,ilutp,droptol,1e-6);
[L,U] ilu(A,options);
[x_precond,fl1,rr1,it1,rv1] bicgstab(A,b,tol,maxit,L,U);
fl1
fl1 0
rr1
rr1 3.8661e-14
it1
it1 3 在第三次迭代中使用 ilu 预条件子产生的相对残差 rr1 小于请求的容差 1e-12。输出 rv1(1) 为 norm(b)输出 rv1(end) 为 norm(b-A*x1)。 可以通过绘制每个二分之一迭代的相对残差来跟踪 bicgstab 的进度。绘制每个解的残差历史记录图并添加一条表示指定容差的线。
semilogy((0:numel(rv0)-1)/2,rv0/norm(b),-o)
hold on
semilogy((0:numel(rv1)-1)/2,rv1/norm(b),-o)
yline(tol,r--);
legend(No preconditioner,ILU preconditioner,Tolerance,LocationEast)
xlabel(Iteration number)
ylabel(Relative residual)
如图所示 提示 此函数返回的不完全分解可用作通过迭代方法例如 bicg、bicgstab 或 gmres求解的线性系统的预条件子。
