回 帖 发 新 帖 刷新版面

主题:新手求助,用幂法求特征值不收敛

请各位大哥指点,小弟首次发帖,这周要交作业又不想抄
我用幂法求一个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

回复列表 (共1个回复)

沙发

我也没做过这个作业.

我来回复

您尚未登录,请登录后再回复。点此登录或注册