24小时热门版块排行榜

返回列表

chaoszero253

新虫 (初入文坛)

应助: 0 (幼儿园)
金币: 5
帖子: 5
在线: 49分钟
虫号: 6900181
注册: 2017-07-04
性别: GG
专业: 光学

[求助] 关于Jacobi法矩阵对角化

想用fortran语言通过Jacobi方法对角化一个对称矩阵，但是求出来的结果和matlab内置的对角化程序员结果不一样。请问我的对角化程序是哪里出问题了。

subroutine SOLVE(A,N,tezheng,tol)
   !jacobi法对角化
   implicit real*8(a-z)
   integer::N
   real*8::A(N,N),tezheng(N)
   real*8::A1(N,N),R(N,N),RT(N,N),U(N,N),thi(N,N)
   integer::i,j,k,m,c,b,z,t
   integer::p,q,x,y
   real*8 sum
   do i=1,N
      do j=1,N
            A1(i,j)=A(i,j)
      end do
   end do
   k=0
   p=0
   m=1
   do while(m==1)
      m=0
      do i=1,N-1
            do j=i+1,N
               R=0
               do q=1,N
                  R(q,q)=1
               end do
               if(abs(A1(i,j))>tol) then
                  m=1
!       判断非对角矩阵元是否大于给出的阈值，如果是则建立旋转矩阵R来归零矩阵元
                  if(A1(i,i)==A1(j,j))then
                        R(i,i)=1/1.41421356
                        R(j,j)=1/1.41421356
                        R(i,j)=1/1.41421356
                        R(j,i)=-1/1.41421356
                  else
                        thi(i,j)=0.5*atan(2*A1(i,j)/(A1(i,i)-A1(j,j)))
                        R(i,i)=cos(thi(i,j))
                        R(j,j)=cos(thi(i,j))
                        R(i,j)=-sin(thi(i,j))
                        R(j,i)=sin(thi(i,j))
                  end if
                  do x=1,N
                        do y=1,N
                           RT(x,y)=R(y,x)
                        end do
                  end do
                  U=matmul(RT,A1)
                  A1=matmul(U,R)
               end if
            end do
      end do
      p=p+1
      print *,p
      if (p==100) exit
   end do
   do i=1,N
      tezheng(i)=A1(i,i)
   end do
end subroutine SOLVE

回复此楼

» 猜你喜欢

深圳大学2026年秋博士招生-物理学-活性胶体方向-高永祥课题组已经有18人回复
论物质与能量的统一模型及物理现象解释已经有0人回复
物理学I论文润色/翻译怎么收费? 已经有101人回复
基于基元I统一理论的数学相关应用推导已经有0人回复
基元I统一理论：宇宙本质、层级演化与修炼文明的本源规律已经有1人回复
基元I理论下三大核心空间现象精准推导与细节解析已经有0人回复
基于基元 I 统一理论的反重力理论推导已经有0人回复
基于基元I统一理论的量子力学本源推导已经有1人回复
推荐一款可以AI辅助写作的Latex编辑器SmartLatexEditor，超级好用，AI润色，全免费已经有20人回复
【EI|Scopus 双检索】第六届智能机器人系统国际会议（ISoIRS 2026）已经有0人回复
2026年第四届电动车与车辆工程国际会议（CEVVE 2026）已经有0人回复

1楼 2018-09-15 10:10:40

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

相关版块跳转我要订阅楼主 chaoszero253 的主题更新

返回列表