左除\和右除/运算符,在上一节中讨论过,使用直接解算器来求解形式为x= A\ b或x= b/ AOctave还包括使用迭代技术求解稀疏线性方程的许多函数。
x = pcg (A, b, tol, maxit, m1, m2, x0, …) ¶x = pcg (A, b, tol, maxit, M, [], x0, …) ¶[x, flag, relres, iter, resvec, eigest] = pcg (A, b, …) ¶求解线性方程组A * x= b采用预处理共轭梯度迭代法。
输入参数为:
Afcn使得Afcn(x)=A*x。的附加参数Afcn可能在之后通过x0.
A必须是埃尔米特和正定(HPD)。如果pcg检测A不是肯定的,打印一个警告旗帜设置输出。
b- A * x。如果标准b- A * x)≤ tol标准b)如果tol则使用1e-6的公差。m= p1 * 第2页使得inv(p1) * A*inv(第2页)是HPD,共轭梯度法正式应用于线性系统inv(p1) * A*inv(第2页) * y=发票(p1) * b具有x=发票(第2页) * y(分裂预处理)。在实践中,在共轭梯度法的每次迭代中,都使系统与矩阵等距m用解决mldivide.如果一个特定的因子分解m= m1 * m2是可用的(例如,的完全Cholesky因式分解a),两个矩阵m1和m2可以通过,并且相对线性系统用mldivide运算符请注意,正确选择预处理器可以显著提高该方法的整体性能。而不是矩阵m1和m2,用户可以传递两个函数,返回应用的逆运算的结果m1和m2到向量。如果m1被省略或为空[],则不应用任何预处理。如果没有的因子分解m是可用的,m2可以省略或保留[],并且输入变量m1可以用来通过预处理器m.
以下参数x0被视为参数,并以适当的方式传递给任何函数(A或m1或m2)已经给予pcg。有关更多详细信息,详见下面的示例。
输出参数为:
A * x= b。如果算法没有收敛,那么x是具有最小残差的迭代。eps标准x2..
长度(resvec)-1.
resvec(我1.是残差的欧几里得范数,并且resvec(我2.是预处理残差形式,在(我-1) 第-次迭代,我= 1, 2, …, iter+1预处理残差范数定义为r * (m\ r)这里的r= b- A * x,详见的说明m如果eigest不是必需的,只是resvec1.返回。eigest1.和最大的eigest2.预处理矩阵的特征值P= m\ A特别地,如果使用norepresenting,则的极端特征值的估计A返回。eigest1.是一个高估和eigest2.是低估了,所以eigest2.eigest1.是的下限cond(P2.,然而在极限中理论上应该等于条件数的实际值。让我们考虑一个三对角矩阵的平凡问题
n=10;A=toeplitz(稀疏([1,1],[1,2],[2,1],1,n));b=A*ones(n,1);M1=ichol(A);#在这种三对角情况下,它对应于chol(A)'M2=M1';M=M1*M2;Afcn=@(x)A*x;Mfcn=@(x)M\x;M1fcn=@(x)M1\x;M2fcn=@(x)M2\x;
示例1:最简单的使用pcg
x=pcg(A,b)
示例2: pcg具有一个计算A * x
x=pcg(Afcn,b)
示例3: pcg具有预处理矩阵M
x=pcg(A,b,1e-06100,M)
示例4: pcg具有作为预处理器的函数
x=pcg(Afcn,b,1e-6100,Mfcn)
示例5: pcg具有预处理矩阵M1和M2
x=pcg(A,b,1e-6100,M1,M2)
示例6: pcg具有作为预处理器的函数
x=pcg(Afcn,b,1e-6100,M1fcn,M2fcn)
示例7: pcg将需要参数的函数作为输入
函数y=Ap(A,x,p)#计算A^p*x y=x;对于i=1:p y=A*y;endfor endfunctionApfcn=@(x,p)Ap(A,x,p);x=pcg(Apfcn,b,[],[],]],[],2],2);
示例8:明确的例子表明pcg使用asplit预处理器
M1=伊科尔(A+0.1*眼(n));#扰动M2=M1'的因子分解;M=M1*M2;##两次迭代后从pcg计算的参考解[x_ref,fl]=pcg(A,b,[],2,M)##分裂预处理[y,fl]=pcg((M1\A)/M2,M1\b,[],2中)x=M2\y#比较x和x_ref
参考文献:
x = pcr (A, b, tol, maxit, m, x0, …) ¶[x, flag, relres, iter, resvec] = pcr (…) ¶求解线性方程组A * x= b利用预处理共轭残差迭代法。
输入参数为
A * x原则上A应该是对称的和非奇异的;如果聚合酶链式反应查找A为了在数字上保持棱角分明,您将get一条警告消息和旗帜将设置outputparameter。b- A * x。如果标准b- A * x) <=tol标准b- A * x0)如果tol为空或被省略,函数集tol=1e-6默认情况下。[]已申请maxit或聚合酶链式反应参数较少,则使用等于20的默认值。聚合酶链式反应 P * x= m\ b具有P= m\ A请注意,正确选择第二种方法可以显著提高该方法的整体性能。而不是矩阵m,用户可以传递一个函数,该函数返回应用的逆的结果m到向量(通常这是使用预处理器的首选方式)。如果[]提供用于m或m省略,不应用修复。以下参数x0被视为参数,并以适当的方式传递给任何函数(A或m)传递给聚合酶链式反应。有关更多详细信息,详见下面的示例。
输出参数为
A * x= b.
旗帜= 0表示收敛的解和下式给出的公差标准tol令人满意。旗帜1.意味着maxit达到了盈利计数的极限。旗帜3.返回a聚合酶链式反应细分,详见[1]。resvec(i)包含残差的欧几里得范数(我-1) 第-次迭代,我= 1,2, …, iter+1.
让我们考虑一个对角矩阵的平凡问题(我们利用a的简约性)
n=10;A=稀疏(diag(1:n));b=兰特(N,1);
示例1:最简单的使用聚合酶链式反应
x=聚合酶链式反应(A,b)
示例2: 聚合酶链式反应具有一个计算A * x.
函数y=apply_a(x)y=[1:10]'.*x;endfunctionx=pcr(“apply_a”,b)
示例3:预处理迭代,具有完整diag函数。预处理器(很奇怪,因为即使是原始矩阵A是平凡的)定义为一个函数
函数y=apply_m(x)k=地板(长度(x)-2);y=x;y(1:k)=x(1:k[1:k]';endfunction[x,flag,relres,iter,resvec]=。。。pcr(A,b,[],[],“apply_m”)半逻辑([1:iter+1],resvec);
示例4:最后,一个依赖于一个参数的预处理器k.
函数y=apply_m(x,varargin)k=varargin{1};y=x;y(1:k)=x(1:k[1:k]';endfunction[x,flag,relres,iter,resvec]=。。。pcr(A,b,[],[],“apply_m”',[],3)参考
W.Hackbusch,大型稀疏方程组的迭代解法,第9.5.4节;施普林格,1994年
使用预处理矩阵可以加快迭代解算器收敛到解的速度M.在这种情况下,线性方程M1.x= M1.A\ b而是被解决。典型的预处理矩阵是原始矩阵的部分因子分解。
L = ichol (A) ¶L = ichol (A, opts) ¶计算稀疏平方矩阵的不完全Cholesky因子分解A.
默认情况下,伊科尔只使用的下三角矩阵A并返回一个较低的三角因子L使得L*L近似A.
该子程序给出的因子可以用作通过迭代方法(如PCG(预处理共轭梯度))求解的线性方程系统的预处理器。
因子分解可以通过传递结构体中的参数来修改opts。参数名称是结构体的一个字段,设置是字段的值。键和值区分大小写。
因子分解的类型。
nofill默认没有填充的不完全Cholesky因子分解(IC(0))。
ict具有阈值下降的不完全Cholesky因子分解(ICT)。
非负标量阿尔法的不完全Cholesky因子分解A+ 阿尔法*diag(diag(A))而不是A。当A不是肯定的。默认值为0。
一种非负标量,指定在执行ICT时因式分解的下降容限。默认值为0,生成完整的Cholesky因子分解。
的非对角线条目L设置为0,除非
abs(L(i,j)>=滴剂*范数(A(j:结束,j),1).
改进的不完全Cholesky因子分解:
关默认行和列的总和不一定保留。
on的对角线L被修改为即使在因子分解过程中删除了元素,也可以保留行(和列)和。保留的关系是:AeL * Le,其中e是一个1的向量。
lower默认只使用的下三角矩阵A并返回一个较低的三角因子L使得L*L近似A.
upper仅使用的上三角矩阵A并返回一个上三角因子U使得U'*U近似A.
示例
下面的问题演示了如何用完全Cholesky分解和完全Choleski分解分解样本对称正定矩阵。
A=[0.37,-0.05,-0.05、-0.07;-0.05,0.116,0.0,-0.05;-0.05、0.0、0.116、-0.05;0.07,-0.05和-0.05,0.202];A=稀疏(A);nnz(tril(A))ans=9L=chol(A,“较低”);nnz(L)ans=10范数(A-L*L',“fro”)/范数(A,“fro”)ans=1.1993e-16opts.type=“nofill”;L=ichol(A,opts);nnz(L)ans=9范数(A-L*L',“fro”)/范数(A,“fro”)ans=0.09736
分解的另一个例子是用于求解单位平方上的边值问题的有限差分矩阵。
nx=400;ny=200;hx=1/(nx+1);hy=1/(ny+1);Dxx=spdiags([ones(nx,1),-2*ones(nx,1;Dyy=spdiags([ones(ny,1),-2*ones(ny+1),ones(nny,1)],[-10 1],ny,ny)/(hy^2);A=-克朗(Dxx,speye(ny))-克朗(speye(nx),Dyy);nnz(tril(A))ans=239400opts.type=“nofill”;L=ichol(A,opts);nnz(tril(A))ans=239400标准(A-L*L',“fro”)/标准(A,“fro”)ans=0.063227
实现算法的参考文献:
[1] Y.萨阿德。“预处理技术。稀疏线性系统的迭代方法,PWS出版公司,1996年。
[2] M.Jones,P.Plassmann:一种改进的不完全Holesky因子分解, 1992.
LUA = ilu (A) ¶LUA = ilu (A, opts) ¶[L, U] = ilu (…) ¶[L, U, P] = ilu (…) ¶计算稀疏平方矩阵的不完全LU因子分解A.
ilu返回单位下三角矩阵L,上三角矩阵U,以及可选的置换矩阵P,使得L*U近似P*A.
该子程序给出的因子可以用作通过迭代方法(如BICG(双共轭梯度)或GMRES(广义最小残差法))求解的线性方程组的渐近性的预处理器。
因子分解可以通过传递结构体中的参数来修改opts。参数名称是结构体的一个字段,设置是字段的值。键和值区分大小写。
类型因子分解的类型。
nofill默认没有填充的ILU因子分解(ILU(0))。
其他支持的参数:密尔.
滑稽克劳特版本的ILU因子分解(ILUC)。
其他支持的参数:密尔, 滴剂.
ilutp具有阈值和枢轴的ILU因子分解。
其他支持的参数:密尔, 滴剂, udiag,脱粒.
滴剂一种非负标量,用于指定因式分解的下降容差。默认值为0,这将返回完整的LU因子分解。
的非对角线条目U设置为0,除非
abs(U(i,j)>=滴剂*范数(Aj.
的非对角线条目L设置为0,除非
abs(L(i,j)>=滴剂*范数(AjU(j,j).
密尔改进的不完全LU因子分解:
一行行和修正的不完全LU因子分解。因子分解保留行和:AeL * Ue,其中e是一个1的向量。
col列和修正的不完全LU因子分解。因子分解保留列和:eAeL * U.
关默认行和列的总和不一定保留。
udiag如果为true,则上三角因下对角线上的任何零都将替换为局部跌落公差droptol*标准(AjU(j,j)。默认值为false。
脱粒因子分解的枢轴阈值。它的范围可以在0(对角线枢轴)和1(默认值)之间,其中列中的最大幅值条目被选择为枢轴。
如果ilu仅用一个输出调用,返回的矩阵为L+ U-speye(尺寸(A))这里的L是单位下三角矩阵和U是上三角矩阵。
在具有两个输出的情况下,ilu返回单位下三角矩阵L和一个上三角矩阵U对于opts类型ilutp,其中一个因素是基于的值进行排列的opts.milu。何时opts.milu==一行, U是列排列的上三角因子。否则L是下排列单元的下三角因子。
如果有三个命名输出,并且opts.milu！=一行,P如此返回L和U是不完全因子P*A什么时候opts.milu==一行, P如此返回L和U是的不完整因素A*P.
示例
A=画廊(“neumann”,1600)+speye(1600);opts.type=“nofill”;nnz(A)ans=7840nnz(lu(A))ans=126478nnz(ilu(A,opts))ans=7840
这表明A具有7840个非零,完全LU因子分解具有126478个非零;不完全LU因子因子分解具有0级填充,具有7840非零,与A。摘自:https://www.mathworks.com/help/matlab/ref/ilu.html
A=画廊(“wathen”,10,10);b=总和(A,2);tol=1e-8;maxit=50;opts.type=“crout”;opts.droptol=1e-4;[L,U]=ilu(A,opts);x=bicg(A,b,tol,maxit,L,U);范数(A*x-b,inf)
本例使用ILU作为随机FEM矩阵的预处理器,该矩阵具有报警条件号。没有L和UBICG不会进行任何转换。
版权所有 © 2024 Octave中文网
ICP备案/许可证号：黑ICP备2024030411号