主题:新手求助,用幂法求特征值不收敛
请各位大哥指点,小弟首次发帖,这周要交作业又不想抄
我用幂法求一个501*501的矩阵,主对角线上元素是ai,两侧是b,再两侧是c
我这样计算为什么不收敛呢?
我抽检了几个u和y的元素,都是趋于零,左乘A应该把向量集中到了最大特征向量的方向啊,怎么算出来的beida越来越大不收敛呢?
program test1
real(kind=8),DIMENSION(1:501,1:501)::A;REAL(kind=8),DIMENSION(1:501)::adj
INTEGER I,J,K
real*8 ::b,c,e,keci
!将A转换为5*501矩阵所需
real*8 ::f2=0,beida=0,beida1=0,wucha
REAL(kind=8),DIMENSION(1:501)::u,y
!幂法迭代使用
PARAMETER (b=0.16,c=-0.064,e=2.718281828459045,keci=10**(-12))
do I=1,501
adj(I)=(1.64-0.024*I)*sin(0.2*I)-0.64*e**(0.1/I)
END do
A(1,1:2)=0;A(2,1)=0;A(4,501)=0;A(5,500:501)=0
A(1,3:501)=c;A(2,2:501)=b;A(3,1:501)=adj(1:501);A(4,1:500)=b;A(5,1:499)=c
!将A转换为5*501矩阵并赋值
u(1:501)=0.02;j=1
do while(wucha>=keci)!当误差不低于keci时进行下一次迭代
j=j+1 !计数器
do i=1,501
f2=f2+u(i)*u(i)
end do
f2=sqrt(f2)
!计算u的2范数
y=u/f2 !把u归一化
do i=3,499
u(i)=adj(i)*y(i)+b*(y(i+1)+y(i-1))+c*(y(i+2)+y(i-2))
end do
u(1)=adj(1)*y(1)+b*y(2)+c*y(3)
u(2)=adj(2)*y(2)+b*(y(1)+y(3))+c*y(4)
u(501)=adj(501)*y(501)+b*y(500)+c*y(499)
u(500)=adj(500)*y(500)+b*(y(499)+y(501))+c*y(498)
!对u左乘A进行迭代
beida1=beida
do i=1,501
beida=beida+y(i)*u(i)
end do !y乘u求beida
wucha=abs(beida-beida1)/abs(beida)!收敛条件
write(*,*) "beida=",beida,"beida1=",beida1,"j=","wucha=",wucha,"j=",j
if(j>2000) write(*,*) "fuck,it is failed",j
if(j>2000) exit
!溢出
end do
end program
我用幂法求一个501*501的矩阵,主对角线上元素是ai,两侧是b,再两侧是c
我这样计算为什么不收敛呢?
我抽检了几个u和y的元素,都是趋于零,左乘A应该把向量集中到了最大特征向量的方向啊,怎么算出来的beida越来越大不收敛呢?
program test1
real(kind=8),DIMENSION(1:501,1:501)::A;REAL(kind=8),DIMENSION(1:501)::adj
INTEGER I,J,K
real*8 ::b,c,e,keci
!将A转换为5*501矩阵所需
real*8 ::f2=0,beida=0,beida1=0,wucha
REAL(kind=8),DIMENSION(1:501)::u,y
!幂法迭代使用
PARAMETER (b=0.16,c=-0.064,e=2.718281828459045,keci=10**(-12))
do I=1,501
adj(I)=(1.64-0.024*I)*sin(0.2*I)-0.64*e**(0.1/I)
END do
A(1,1:2)=0;A(2,1)=0;A(4,501)=0;A(5,500:501)=0
A(1,3:501)=c;A(2,2:501)=b;A(3,1:501)=adj(1:501);A(4,1:500)=b;A(5,1:499)=c
!将A转换为5*501矩阵并赋值
u(1:501)=0.02;j=1
do while(wucha>=keci)!当误差不低于keci时进行下一次迭代
j=j+1 !计数器
do i=1,501
f2=f2+u(i)*u(i)
end do
f2=sqrt(f2)
!计算u的2范数
y=u/f2 !把u归一化
do i=3,499
u(i)=adj(i)*y(i)+b*(y(i+1)+y(i-1))+c*(y(i+2)+y(i-2))
end do
u(1)=adj(1)*y(1)+b*y(2)+c*y(3)
u(2)=adj(2)*y(2)+b*(y(1)+y(3))+c*y(4)
u(501)=adj(501)*y(501)+b*y(500)+c*y(499)
u(500)=adj(500)*y(500)+b*(y(499)+y(501))+c*y(498)
!对u左乘A进行迭代
beida1=beida
do i=1,501
beida=beida+y(i)*u(i)
end do !y乘u求beida
wucha=abs(beida-beida1)/abs(beida)!收敛条件
write(*,*) "beida=",beida,"beida1=",beida1,"j=","wucha=",wucha,"j=",j
if(j>2000) write(*,*) "fuck,it is failed",j
if(j>2000) exit
!溢出
end do
end program