numerical_analysis1.m

CG方法 高等数值分析 PCG方法 最速下降法

clear;
n=1000;
k=400;
H=hilb(n);
%Get A5----------------------------------------
d0=diag(H);
d1=diag(H,-1);
d2=diag(H,1);
d3=diag(H,-2);
d4=diag(H,2);
A3=diag(d0)+diag(d1,-1)+diag(d2,1);
A5=diag(d0)+diag(d1,-1)+diag(d2,1)+diag(d3,-2)+diag(d4,2);

xs=10.*rand(n,1);%[6.0854 0.1576 0.1635 1.9007 5.8692 0.5758 3.6757 6.3145 7.1763 6.9267 0.8408 4.5436 4.4183 3.5325].';%10.*rand(n,1);

b=H*xs;
tic
%PCG Algorithms*********************************
%initialization--------------------------------
R=zeros(n,k+1);
X=zeros(n,k+1);
P=zeros(n,k+1);
Z=zeros(n,k+1);
alfa=zeros(1,k);
beta=zeros(1,k);
r=zeros(1,k+1);
R(:,1)=b;
P(:,1)=A5\b;
Z(:,1)=P(:,1);
r(1)=normest(b);
%iteration-------------------------------------
for j=1:k
    alfa(j)=( R(:,j).'*Z(:,j) )./( P(:,j).'*(H*P(:,j)) );
    X(:,j+1)=X(:,j)+alfa(j).*P(:,j);
    R(:,j+1)=R(:,j)-alfa(j).*(H*P(:,j));
    r(j+1)=normest(R(:,j+1));
    if r(j+1)<=1e-7
        break;
    end
    Z(:,j+1)=A5\R(:,j+1);
    beta(j)=( R(:,j+1).'*Z(:,j+1) )./( R(:,j).'*Z(:,j) );
    P(:,j+1)=Z(:,j+1)+beta(j).*P(:,j);
end
x5=X(:,j+1);
step5=j
toc
tic
%PCG Algorithms*********************************
%initialization--------------------------------
R=zeros(n,k+1);
X3=zeros(n,k+1);
P=zeros(n,k+1);
Z=zeros(n,k+1);
alfa=zeros(1,k);
beta=zeros(1,k);
r=zeros(1,k+1);
R(:,1)=b;
P(:,1)=A3\b;
Z(:,1)=P(:,1);
r(1)=normest(b);
%iteration-------------------------------------
for j=1:k
    alfa(j)=( R(:,j).'*Z(:,j) )./( P(:,j).'*(H*P(:,j)) );
    X3(:,j+1)=X3(:,j)+alfa(j).*P(:,j);
    R(:,j+1)=R(:,j)-alfa(j).*(H*P(:,j));
    r(j+1)=normest(R(:,j+1));
    if r(j+1)<=1e-7
        break;
    end
    Z(:,j+1)=A3\R(:,j+1);
    beta(j)=( R(:,j+1).'*Z(:,j+1) )./( R(:,j).'*Z(:,j) );
    P(:,j+1)=Z(:,j+1)+beta(j).*P(:,j);
end
x3=X3(:,j+1);
step3=j
toc
del3=delta(x3,xs)
del5=delta(x5,xs)