高等数值分析作业-第二次实验

高等数值分析第二次实验作业

注：矩阵阶数均为1000阶

T1. 构造例子特征值全部在右半平面时，观察基本的Arnoldi方法和GMRES方法的数值性态，和相应重新启动算法的收敛性。

Answer：

➢关于特征值均在右半平面的矩阵A：

首先构造对角矩阵D，其中D矩阵的对角是由以下形式的2×2的子矩阵组成其对角元：

S=[a−b

b a

]，其中a，b为实数

此时，S矩阵的特征值分别为a+bi和a-bi，这样D=diag(S1,S2,S3……S n)矩阵的特征值均分布在右半平面。另矩阵A=Q T AQ，则A矩阵的特征值也均在右半平面，生成矩阵的过程代码如下所示：

N=500 %生成的矩阵为2N阶

A=zeros(2*N);

%delta控制特征值分布的密集程度，即控制条件数大小，delta越大条件数越大

delta=0.1;

%构造D矩阵

for j=1:N

A(2*j-1,2*j-1)=N+j*delta;

A(2*j-1,2*j)=-N-j*delta;

A(2*j,2*j-1)=N+j*delta;

A(2*j,2*j)=N+j*delta;

end

U = orth(rand(2*N,2*N));

A = U'*A*U;

➢首先进行观察基本的Arnoldi方法和GMRES方法的数值性态，考虑N=500，即矩阵为1000阶，取X0=zeros(N,1)，b=ones(N,1)，收敛准则为e=10-6;

Arnoldi方法函数如下：

function [xm,error,num]=Arnoldi(A,x0,b,e) N=size(A,1);

r=b-A*x0;

belta=norm(r);

v=r/belta;

V=v;

j=0;

while norm(r)>e & j

j=j+1;

num(j)=j;

w=A*V(1:N,j);

for i=1:j h(i,j)=V(1:N,i)'*w; w=w-h(i,j)*V(1:N,i);

end

h(j+1,j)=norm(w);

if h(j+1,j) ~=0

v=w/h(j+1,j);

V=[V v];

end

e1=zeros(j,1);

e1(1)=1;

be1=belta*e1;

try

[L,U,S]=lu(h(1:j,1:j)); be1=S*be1;

lym=LX(L,be1);

ym=UX(U,lym);

xm=x0+V(1:N,1:j)*ym; r=b-A*xm;

end

error(j)=log10(norm(r)); end

end

GMRES方法的函数如下：

function [xm,error,num] = GMRES(A,x0,b,e) %ARNOLDI Summary of this function goes here

% Detailed explanation goes here

N=size(A,1);

r=b-A*x0;

belta=norm(r);

v=r/belta;

V=v;

j=0;

er=1000;

while er > e & j

j=j+1;

num(j)=j;

w=A*V(1:N,j);

for i=1:j

h(i,j)=V(1:N,i)'*w; w=w-h(i,j)*V(1:N,i);

end

h(j+1,j)=norm(w);

if h(j+1,j) ~=0

v=w/h(j+1,j);

V=[V v];

end

try

[Q,R]=qr(h(1:j+1,1:j)); gk=Q'*belta*eye(j+1,1); ym=minresYK(R,gk);

xm=x0+V(1:N,1:j)*ym;

er=gk(j+1);

end

error(j)=log10(er);

end

end

其中UX，LX和minresYK是自己编写的求解（上、下）三角矩阵的函数，代码如下：

function yk=LX(TK,b) n=size(TK,1);

yk=zeros(n,1);

yk(1)=b(1)/TK(1,1); for i=2:n

yk(i)=b(i);

for j=1:i-1

yk(i)=yk(i)-TK(i,j)*yk(j);

end

yk(i)=yk(i)/TK(i,i); end

end function yk=UX(TK,b)

n=size(TK,2);

yk=zeros(n,1);

yk(n)=b(n)/TK(n,n);

for i=1:n-1

k=n-i;

yk(k)=b(k);

for j=k+1:n

yk(k)=yk(k)-

TK(k,j)*yk(j);

end

yk(k)=yk(k)/TK(k,k);

end

end

function yk=minresYK(TK,b)

n=size(TK,2);

yk=zeros(n,1);

yk(n)=b(n)/TK(n,n);

for i=1:n-1

k=n-i;

yk(k)=b(k);

for j=k+1:n

yk(k)=yk(k)-

TK(k,j)*yk(j);

end

yk(k)=yk(k)/TK(k,k);

end

end

两种方法的计算结果如下所示：