FOM-GMRES-CG等Krylov子空间迭代算法以及Matlab程序

自我o救赎

已于 2023-04-20 11:01:17 修改

阅读量1.3k

点赞数 3

分类专栏：数学大类专栏文章标签：算法 matlab Krylov方法

于 2023-04-19 08:27:03 首次发布

本文链接：https://blog.csdn.net/y15520833229/article/details/130235246

版权

数学大类专栏专栏收录该内容

146 篇文章 8 订阅

订阅专栏

在这里插入图片描述

算法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

在这里插入图片描述

算法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

在这里插入图片描述
算法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循环语句，只有当j>=2时，才会执行
    for i=1:j-1
        H(i:i+1,j)=[c(i),s(i);-s(i),c(i)]*H(i:i+1,j); 
    end
    %因此前面这两段程序在原文的算法中，顺序相反，但并无影响
    %见我的博客《Krylov子空间迭代算法》  https://blog.csdn.net/y15520833229/article/details/130165391
    if H(j+1,j)==0
        m=j;
        break;
    else
        V(:,j+1)=w(:,j)/H(j+1,j);
    end
     %向量c和s的计算赋值当j=1，就从这里开始
    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

%%%调用如下例子进行测试
>>A=[ 20, -10.,   0.,   0.,   0.,   0.,   0.,   0.,   0.;
 -10.,  20., -10.,   0.,   0.,   0.,   0.,   0.,   0.;
  0., -10.,  20., -10.,   0.,   0.,   0.,   0.,   0.;
  0.,   0., -10.,  20., -10.,   0.,   0.,   0.,   0.;
 0. ,  0.,   0., -10.,  20., -10.,   0.,   0.,   0.;
  0.,   0.,   0.,   0., -10.,  20., -10.,   0.,   0.;
 0.,   0.,   0.,   0.,   0., -10.,  20., -10.,   0.;
  0.,   0.,   0.,   0.,   0.,   0., -10.,  20, -10.;
  0.,   0.,   0.,   0.,   0.,   0.,   0., -10.,  20.];
>>b=[0.11407755, 0.048673,   0.03304509, 0.02532112, 0.02064553, 0.01749003, 0.01520812, 0.01347669, 0.0121155 ]';
>>x0=ones(9,1);
>> GMRES_economic(A,x0,2*10^(-8),b)

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

在这里插入图片描述
还见博客
https://blog.csdn.net/xiexinning/article/details/121880965?app_version=5.15.0&code=app_1562916241&csdn_share_tail=%7B%22type%22%3A%22blog%22%2C%22rType%22%3A%22article%22%2C%22rId%22%3A%22121880965%22%2C%22source%22%3A%22y15520833229%22%7D&uLinkId=usr1mkqgl919blen&utm_source=app