JACOBI-EIGEN

 ! Program name:JACOBI-EIGEN

real, allocatable :: A(:,:), O(:,:),OT(:,:), M(:,:)

real :: alfa

integer :: i, j, k, l

 

Print *, "Enter the dimension of the symmetric matrix A"

Read *, n

 

Allocate(A(n,n))

Allocate(O(n,n))

Allocate(OT(n,n))

Allocate(M(n,n))

 

Do i=1,n

   Do j=1,n

     Print *, 'A(',i,',',j,')= '

     Read *, A(i,j)

   End do

end do

do i=1,n

     write (*,*) (A(i,j),j=1,n)

end do

print *,

 

Do i=1,n

   Do j=1,n

   If (i .eq. j) O(i,j)=1

   if (i .ne. j) O(i,j)=0

   end do

end do

OT=transpose(O)

 

do i=1,n-1

do j=i+1,n

 

       alfa=0.5*atan(2*A(i,j)/(A(i,i)-A(j,j)))

 

       O(i,i)=cos(alfa)

       O(j,j)=cos(alfa)

       O(i,j)=-sin(alfa)

       O(j,i)=sin(alfa)

       OT(i,i)=cos(alfa)

       OT(j,j)=cos(alfa)

       OT(i,j)=sin(alfa)

       OT(j,i)=-sin(alfa)

 

      M=matmul(A,O)

      A=matmul(OT,M)

 

      Do k=1,n

        Do l=1,n

         If (k .eq. l) O(k,l)=1

         if (k .ne. l) O(k,l)=0

        end do

      end do

      OT=transpose(O)

 

    end do

end do

 

do i=1,n

         write (*,*) (A(i,j),j=1,n)

end do

print *

 

stop

end

Comentários

Postar um comentário

Postagens mais visitadas deste blog

HYDROGEN-RADIAL

!--------------------------------------------------------------- !Program name: HYDROGEN-RADIAL !Radial equation of the hydrogen-like atom !Objective: plot the radial wave function for each n,l pair ! The equation is atomic units !--------------------------------------------------------------- !The Coulomb potential energy equation is: ! V(r) = -2Z/r ! The ceentrifugal potential equation is: ! VL=l(l+1)/r^2 ! The energy is: ! E=-2*Z/n^2 !--------------------------------------------------------------- ! Use of Numerov method to find the wave function ! Function fn from Numerov equation is ! Fn = 1 + K2*h^2/12 , where K2=-2(Veff - E) ! In the code, h is changed into dr !---------------------------------------------------------------   Program hydrogen   Implicit none   integer :: mesh, n, l, i, j, zeta, m, maxiter=100   Real*8 rmin, dr, zmesh, y_out_m   double precision :: e, fh12, norm, eps=1.0E-6    double precision, allocatable :: r(:), r2
Leia mais

HYDROGEN-POTENTIAL

  !--------------------------------------------------------------- !Program name: HYDROGEN-POTENTIAL !Objective: plot the effective potential energy as a function of radial distance !--------------------------------------------------------------- !The Coulomb potential energy equation is: ! V(r) = -KZq 2 /r = -(Kq 2 )(Z/r) ! K(Coulomb constant)= 9*10 9 Nm 2 C -2 ! q = - 1.602 x 10 -19 C ! Kq 2 =2.3 x 10 -28 ! The unit is Joule. We have to change into Rydberg. !---------------------------------------------------- ! The V(centrifugal) equation is: ! Vcent=6.1 x 10 -39 l(l+1)/r 2 (see text) ! The unit is Joule. We have to change into Rydberg !--------------------------------------------------------------- ! The unit of the radial distance, r, is meter. !---------------------------------------------------------------   Program hydrogen Implicit none   integer :: mesh, n, l, i, zeta   Real*8 rmin, dr, rmax   double precision :: oldr, e, kq2, h2m
Leia mais

ONE_PART_QHO

  !--------------------------------------------------------------- !Program name: ONE_PART_QHO !Original code: Giannozzi,P., Ercolessi, F. and Gironcoli, S. !Modified by Caio Lima Firme !One particle quantum harmonic oscillator !PARABOLIC POTENTIAL WELL !Objective: to plot the wave function for each node !--------------------------------------------------------------- !The dimensionless physical quantities: !Potential energy: V(x)=0.5X^2    ! Energy: E=1/2+n !--------------------------------------------------------------- ! Use of Numerov method to find the wave function ! Function fn from Numerov equation is ! Fn = 1 + K2*h^2/12 , where K2=2m/((H/2pi)^2)*(E-V) ! H=Planck constant !-------------------------------------------------------------- -   integer :: grid, i, icl   integer :: n, hn, ncross, j, iter   double precision :: xmax, h, fh12, norm, djump, fac   double precision :: Vupper, Vlower, e, k   double precision :: x(0:300), y(0:300), V
Leia mais