文章目录
krylov子空间代码
投影方法
给定初值 x ( 0 ) ∈ R n x^{(0)}\in \mathbb{R}^n x(0)∈Rn,采用仿射空间 x ( 0 ) + K x^{(0)}+\mathcal{K} x(0)+K,其中 K \mathcal{K} K 是krylov子空间。找到一个 x ~ ∈ x ( 0 ) + K \tilde{x}\in x^{(0)}+\mathcal{K} x~∈x(0)+K,令 x ~ = x ( 0 ) + x ^ \tilde{x}=x^{(0)}+\hat{x} x~=x(0)+x^ 使得 r 0 − A x ^ ⊥ L (1) r_0-A\hat{x}\perp\mathcal{L}\tag{1} r0−Ax^⊥L(1)此时可以证明 x ~ \tilde{x} x~ 是 $ Ax=b$ 的近似解。
krylov子空间
A ∈ R n × n , r ∈ R n , K m ( A , r ) ≜ s p a n { r , A r , . . . , A m − 1 r } ⊆ R n A\in \mathbb{R}^{n\times n},r\in\mathbb{R}^n,\mathcal{K}_m(A,r)\triangleq span\{r,Ar,...,A^{m-1}r\}\subseteq\mathbb{R}^n A∈Rn×n,r∈Rn,Km(A,r)≜span{r,Ar,...,Am−1r}⊆Rn
求krylov子空间的方法:Arnoldi过程
Arnoldi过程是一种特殊的求正交基的方法,求出来的基是 V m = s p a n { r , A r , . . . , A ( m − 1 ) r } V_m=span\{r,Ar,...,A^{(m-1)}r\} Vm=span{r,Ar,...,A(m−1)r},同时会计算出 H m + 1 , m H_{m+1,m} Hm+1,m的矩阵, H m = H ( 1 : m , m ) H_m=H(1:m,m) Hm=H(1:m,m)是Hessenberg矩阵。 V m V_m Vm和 H H H之间的关系如下: A V m = V m + 1 H m + 1 , m = V m H m , m + h m + 1 , m v m + 1 e m T AV_m=V_{m+1}H_{m+1,m}=V_{m}H_{m,m}+h_{m+1,m}v_{m+1}e_m^T AVm=Vm+1Hm+1,m=VmHm,m+hm+1,mvm+1emT V m T A V m = H m , V m T r = V m T β v 1 = β e 1 , β = ∥ r ∥ 2 (2) V_m^TAV_m=H_m, V_m^Tr=V_m^T\beta v_1=\beta e_1, \beta=\Vert{r}\Vert_2\tag{2} VmTAVm=Hm,VmTr=VmTβv1=βe1,β=∥r∥2(2)
算法1:基于Gram-Schmidt正交化的Arnoldi过程
function [V,H]=arnoldi_GS(A,r,m)
%输入初始的向量r,K维数m,原始矩阵A;
%输出一组标准正交基Vm,hessenberg matrix H;
%一种特殊的求正交基的方法,求出来的基是span{r,Ar,...,A^(m-1)r};
n=size(A,1);
V=zeros(n,m+1);
H=sparse(m+1,m);
V(:,1)=r/norm(r);
for j=1:m
w=A*V(:,j);
for i=1:j
H(i,j)=dot(w,V(:,i));
end
for i=1:j
w=w-H(i,j)*V(:,i);
end
H(j+1,j)=norm(w);
if H(j+1,j)==0
break;
else
V(:,j+1)=w/H(j+1,j);
end
end
end
算法2:基于MGS的Arnoldi过程
function [V,H]=arnoldi_GMS(A,r,m)
%输入初始的r,和K维数m,原始矩阵A;
%输出一组标准正交基Vm,hessenberg matrix H;
%一种特殊的求正交基的方法,求出来的基是span{r,Ar,...,A^(m-1)r};
n=size(A,1);
V=zeros(n,m+1);
H=sparse(m+1,m);
V(:,1)=r/norm(r);
for j=1:m
w=A*V(:,j);
for i=1:j
H(i,j)=dot(w,V(:,i));
w=w-H(i,j)*V(:,i);
end
H(j+1,j)=norm(w);
if H(j+1,j)==0
break;
else
V(:,j+1)=w/H(j+1,j);
end
end
end
计算近似解
设 K m = { v 1 , . . . , v m } , L m = { w 1 , . . . , w m } , V = [ v 1 , . . . , v m ] , W = [ w 1 , . . . , w m ] \mathcal{K}_m=\{v_1,...,v_m\},\mathcal{L}_m=\{w_1,...,w_m\},V=[v_1,...,v_m],W=[w_1,...,w_m] Km={v1,...,vm},Lm={w1,...,wm},V=[v1,...,vm],W=[w1,...,wm]
x ^ = V y , x ~ = x ( 0 ) + x ^ (3) \hat{x}=Vy,\tilde{x}=x^{(0)}+\hat{x}\tag{3} x^=Vy,x~=x(0)+x^(3)由(1)知 r 0 − A V y ⊥ w i , i = 1 , 2... , m r_0-AVy\perp w_i,i=1,2...,m r0−AVy⊥wi,i=1,2...,m写成矩阵形式 W T A V y = W T r 0 (4) W^TAVy=W^Tr_0\tag{4} WTAVy=WTr0(4)解出 y y y 带入(2)就是近似解。
FOM方法
我们令 L m = K m , W = V \mathcal{L}_m=\mathcal{K}_m,W=V Lm=Km,W=V,(4)式变成 V T A V y = V T r 0 V^TAVy=V^Tr_0 VTAVy=VTr0, 由(2)知 y = β H m − 1 e 1 y=\beta H_m^{-1}e_1 y=βHm−1e1,故 x ~ = x ( 0 ) + β V m H m − 1 e 1 \tilde{x}=x^{(0)}+\beta V_m H_m^{-1}e_1 x~=x(0)+βVmHm−1e1
残差估计: ∥ r ~ ∥ 2 = h m + 1 , m ∣ y ( m ) ∣ \Vert{\tilde{r}}\Vert_2=h_{m+1,m}|y(m)| ∥r~∥2=hm+1,m∣y(m)∣ 可以利用残差在循环中判断是否达到要求的解,这种方法是实用FOM
算法3:FOM
function [x,error]=FOM(A,x_0,m,b)
%FOM方法 初值和较小的m后就能算出来
%输出近似解x(n维)
%输入矩阵A,b,初值r,K空间的维数m;
r=b-A*x_0;
beta=norm(r);
[V,H]=arnoldi_GS(A,r,m);
e1=eye(m,1);
y=H(1:m,:)\(beta*e1);
error=H(m+1,m)*abs(y(m,1));
x=x_0+V(:,1:size(y))*y;
end
算法4:实用FOM
function [x,error,m]=FOM_economic(A,x_0,e,b)
%FOM方法 给定误差 就能算出来
%输出近似解x(n维)
%输入矩阵A,b,初值r,误差e,K空间的初始维数m;
m=0;
r=b-A*x_0;
beta=norm(r);
V(:,1)=r/norm(r);
while(1)
m=m+1;
w=A*V(:,m);
for i=1:m
H(i,m)=dot(w,V(:,i));
w=w-H(i,m)*V(:,i);
end
H(m+1,m)=norm(w);
e1=eye(m,1);
y=H(1:m,:)\(beta*e1);
error=H(m+1,m)*abs(y(m));
if error<e
break;
else
V(:,m+1)=w/H(m+1,m);
end
end
x=x_0+V(:,1:size(y))*y;
end
GMRES方法
我们令 L m = A K m , W = A V \mathcal{L}_m=A\mathcal{K}_m,W=AV Lm=AKm,W=AV,(4)式变成当且仅当极小化问题的解 m i n x ∈ x ( 0 ) ∥ b − A x ∥ 2 min_{x\in x^{(0)}}\Vert{b-Ax}\Vert_2 minx∈x(0)∥b−Ax∥2
算法5:实用GMRES
function [x,error,m]=GMRES_economic(A,x_0,e,b)
%给定误差 就能算出来
%输出近似解x(n维)
%输入矩阵A,b,初值r,误差e,K空间的初始维数m;
j=0;
r=b-A*x_0;
beta=norm(r);
V(:,1)=r/norm(r);
Eps=beta*eye(size(x_0,1),1);
while(1)
j=j+1;
w(:,j)=A*V(:,j);
for i=1:j
H(i,j)=dot(w(:,j),V(:,i));
w(:,j)=w(:,j)-H(i,j)*V(:,i);
end
H(j+1,j)=norm(w(:,j));
%以上是arnoldi过程
%对H做givens变换 第几行就要做几次givens变换,第一次是给前两行, 最后givens变换是看最新生成的H(j,j)
for i=1:j-1
H(i:i+1,j)=[c(i),s(i);-s(i),c(i)]*H(i:i+1,j);
end
if H(j+1,j)==0
m=j;
break;
else
V(:,j+1)=w(:,j)/H(j+1,j);
end
if abs(H(j,j))>abs(H(j+1,j))
tau=H(j+1,j)/H(j,j);
c(j)=1/sqrt(1+tau^2);
s(j)=c(j)*tau;
else
tau=H(j,j)/H(j+1,j);
s(j)=1/sqrt(1+tau^2);
c(j)=s(j)*tau;
end
H(j,j)=c(j)*H(j,j)+s(j)*H(j+1,j);
H(j+1,j)=0;
Eps(j:j+1)=[c(j),s(j);-s(j),c(j)]*[Eps(j);0];
error=abs(Eps(j+1))*beta;
if abs(Eps(j+1))*beta<e
m=j;
break
end
end
y=H(1:m,1:m)\(Eps(1:m));
x=x_0+V(:,1:size(y))*y;
end
算法6:共轭梯度法
function [x,m, r2] = CG(A, b, error, max_n)
n = size(A, 1);
x = zeros(n, 1);
r = b - A * x;
p = r;
r2 = r' * r;
r2_eps = r2 * error * error;
%r2_old = r2;
m = 1;
r_norm(m) = norm(r);
while ((m < max_n) && (r2 > r2_eps))
s = A * p;
alpha = r2 / (p' * s);
x = x + alpha * p;
r = r - alpha * s;
r2_old = r2;
r2 = r' * r;
beta = r2 / r2_old;
p = r + beta * p;
m = m + 1;
r_norm(m) = norm(r, 2);
end
end
算例结果
testA的规模是
2304
×
2304
2304\times2304
2304×2304 ,当误差小于
1
e
−
9
1e-9
1e−9结束迭代