回 帖 发 新 帖 刷新版面

主题:[讨论]关于数据统计

关于数据统计
现在有一组数据,我想统计下在两个数字间隔内有多少个数据。数据如下:

eigenvalue=  -2.59461975097656     
 eigenvalue=  -2.59461951255798     
 eigenvalue=  -2.59461903572083     
 eigenvalue=  -2.59461879730225     
 eigenvalue=  -2.58745384216309     
 eigenvalue=  -2.58745336532593     
 eigenvalue=  -2.58745312690735     
 eigenvalue=  -2.58745288848877     
 eigenvalue=  -1.87168383598328     
 eigenvalue=  -1.87168371677399     
 eigenvalue=  -1.87168371677399     
 eigenvalue=  -1.87168347835541     
 eigenvalue=  -1.80774712562561     
 eigenvalue=  -1.80774688720703     
 eigenvalue=  -1.80774688720703     
 eigenvalue=  -1.80774676799774     
 eigenvalue=  -1.41421377658844     
 eigenvalue=  -1.41421353816986     
 eigenvalue=  -1.41421353816986     
 eigenvalue=  -1.41421341896057     
 eigenvalue=  -1.30447649955750     
 eigenvalue=  -1.30447638034821     
 eigenvalue=  -1.30447638034821     
 eigenvalue=  -1.30447626113892     
 eigenvalue= -0.316583961248398     
 eigenvalue= -0.316583901643753     
 eigenvalue= -0.316583812236786     
 eigenvalue= -0.316583722829819     
 eigenvalue= -2.427690048989462E-007
 eigenvalue= -1.448813407023408E-007
 eigenvalue= -4.662113894937647E-008
 eigenvalue= -2.762132211842072E-008
 eigenvalue=  2.762132034206388E-008
 eigenvalue=  4.662114605480383E-008
 eigenvalue=  1.448813691240503E-007
 eigenvalue=  2.427690048989462E-007
 eigenvalue=  0.316583663225174     
 eigenvalue=  0.316583782434464     
 eigenvalue=  0.316583901643753     
 eigenvalue=  0.316583961248398     
 eigenvalue=   1.30447638034821     
 eigenvalue=   1.30447649955750     
 eigenvalue=   1.30447649955750     
 eigenvalue=   1.30447673797607     
 eigenvalue=   1.41421353816986     
 eigenvalue=   1.41421365737915     
 eigenvalue=   1.41421377658844     
 eigenvalue=   1.41421389579773     
 eigenvalue=   1.80774664878845     
 eigenvalue=   1.80774712562561     
 eigenvalue=   1.80774712562561     
 eigenvalue=   1.80774724483490     
 eigenvalue=   1.87168335914612     
 eigenvalue=   1.87168371677399     
 eigenvalue=   1.87168383598328     
 eigenvalue=   1.87168395519257     
 eigenvalue=   2.58745336532593     
 eigenvalue=   2.58745360374451     
 eigenvalue=   2.58745384216309     
 eigenvalue=   2.58745455741882     
 eigenvalue=   2.59461951255798     
 eigenvalue=   2.59461951255798     
 eigenvalue=   2.59461975097656     
 eigenvalue=   2.59462022781372     

给定一个参数Ef,让它从-2.6变化到2.6,每隔0.01变化一次,看下再每个间隔内有多少数据。比如在-2.6到-2.59之间,数据数为0;在-2.59到-2.58之间,数据数为4;在-2.58到-2.57之间数据数为4,等等。
我是这样做的,但是没能通过 程序。定义一个数组ds存储数据数,一个数组Ef存储数据间隔。


integer,dimension(:)::ds(1000)  !!ds定义的大些
real,dimension(:)::Ef(1000)  !Ef定义的也大些,便于存放较多数据

ds=0
 do j=1,4*m*n
  do i=1,520
    Ef(i)=2.6-0.01*i
   if(Ef(i)<=eigenvalue(j)<Ef(i+1)) then
   ds(i)=ds(i)+1
    
   end if
   
   open(10,file='taimidu.txt')
     write(10,*) Ef ds(i)
   close(10)
   end do
 end do

再把程序贴上来吧:
program main
  use IMSL
!  use linear_operators
  implicit none
  real, parameter::pi=3.14159265
  complex, parameter:: fi=(0.0, 1.0)
  real t,f,w,temp

  integer,dimension(:)::ds(1000)
  real,dimension(:)::Ef(100)
  complex*8,allocatable::h(:,:)
  complex*8,allocatable::h0(:,:)
  complex*8,allocatable::r0(:,:)
  complex*8,allocatable::l0(:,:)
  
  complex*8,allocatable::eigenvector(:,:)


  real*8, allocatable:: eigenvalue(:)  !特征值的矩阵
  integer ::i,n,m,j,p
  open(1,file='hmat.txt')
  n=4
  m=4 
  
  allocate(h(4*m*n,4*m*n))
  allocate(h0(4*n,4*n),r0(4*n,4*n),l0(4*n,4*n))
  allocate(eigenvalue(4*n*m)) 
  allocate(eigenvector(4*n*m,4*m*n))

  w=0.0
  t=1.0
  f=0.25

  call hmatr(h, n,m)  
 do  i=1,4*n*m
   do j=1,4*n*m

  write(1,*) h(j,i)

   end do
 end do
close(1)


  eigenvalue=eig(h,w=eigenvector) 
  
  open(10,file='eeig.txt')
  do j=1,4*n*m-1
     p=j
    do i=j+1,4*m*n
       if(eigenvalue(i)<eigenvalue(p)) p=i
    end do
     temp=eigenvalue(j)
     eigenvalue(j)=eigenvalue(p)
     eigenvalue(p)=temp
  end do
  
  do i=1,4*m*n
 
   write(10,*) 'eigenvalue=',eigenvalue(i)
   end do
  close(10)
ds=0
 do j=1,4*m*n
  do i=1,51
    Ef(i)=2.6-0.1*i
   if(Ef(i)<=eigenvalue(j)<Ef(i+1)) then
   ds(i)=ds(i)+1
    
   end if
   
   open(10,file='taimidu.txt')
     write(10,*) Ef ds(i)
   close(10)
   end do
 end do
 

 stop
 contains
   subroutine hmatr(h,n,m) 
   integer n,m
   complex h(4*m*n,4*m*n)
   complex h0(4*n,4*n)
   complex r0(4*n,4*n)
   complex l0(4*n,4*n)

   integer:: j,n0,x
   complex:: c(4,4),a(4,4),b(4,4),r(4,4),l(4,4)
 
   h0=0.0
   a=0.0
   b=0.0
   h=0.0
   r0=0.0
   l0=0.0
   a(1, 4)=t
   b=transpose(a)
   c=0.0
   r=0.0
   l=0.0

   
      n0=1
   do while(n0<=2*n-1)

      r(2,1)=t*exp(fi*(pi-2.0*pi*f/3.0*(3.0*n0/2.0+0.25)))
      r(3,4)=t*exp(fi*(pi-2.0*pi*f/3.0*(3.0*(n0+1)/2.0+0.25)))
      l(1,2)=t*exp(-fi*(pi-2.0*pi*f/3.0*(3.0*n0/2.0+0.25)))!conjg(r0(2,1))
      l(4,3)=t*exp(-fi*(pi-2.0*pi*f/3.0*(3.0*(n0+1)/2.0+0.25)))!conjg(r0(3,4))


      r0(2*n0,2*n0-1)=r(2,1)
      r0(2*n0+1,2*n0+2)=r(3,4)
      l0(2*n0-1,2*n0)=l(1,2)
      l0(2*n0+2,2*n0+1)=l(4,3)


      c(1,1)=w
      c(2, 1) = t*exp(fi*(-pi+2.0*pi/3.0*f*(1.5*n0+0.25))) 
      c(1, 2) = t*exp(-fi*(-pi+2.0*pi/3.0*f*(1.5*n0+0.25)))
      
      c(2, 2) = w
      c(2, 3) = t
      c(3, 2) = t
      c(3, 3) = w
      c(3, 4) = t*exp(fi*(-pi+2.0*pi/3.0*f*(1.5*(n0+1)+0.25)))
      c(4, 3) = t*exp(-fi*(-pi+2.0*pi/3.0*f*(1.5*(n0+1)+0.25)))
      c(4, 4) = w
    

       h0(2*n0-1,2*n0-1)=c(1,1)
       h0(2*n0,2*n0-1)=c(2,1)
       h0(2*n0-1,2*n0)=c(1,2)
       h0(2*n0,2*n0)=c(2,2)
       h0(2*n0,2*n0+1)=c(2,3)
       h0(2*n0+1,2*n0)=c(3,2)
       h0(2*n0+1,2*n0+1)=c(3,3)
       h0(2*n0+1,2*n0+2)=c(3,4)
       h0(2*n0+2,2*n0+1)=c(4,3)
       h0(2*n0+2,2*n0+2)=c(4,4)   

       if(2*n0-1==1) then

        h0(2*n0-1:2*n0+2,2*n0+3:2*n0+6)=b
        h0(2*n0-1:2*n0+2,4*n-3:4*n)=a

         else if(2*n0-1<4*n-3)  then

          h0(2*n0-1:2*n0+2,2*n0+3:2*n0+6)=b
          h0(2*n0-1:2*n0+2,2*n0-5:2*n0-2)=a

         else 

            h0(4*n-3:4*n,1:4)=b
            h0(2*n0-1:2*n0+2,2*n0-5:2*n0-2)=a
       end if

       n0=n0+2


   end do

 do x=1,4*(m-1)*n+1,4*n
     h(x:x+4*n-1,x:x+4*n-1)=h0
        if( x==1) then
          h(x:x+4*n-1,x+4*n:8*n)=r0
          h(x:x+4*n-1,(m-1)*4*n+1:m*4*n)=l0
        elseif( x>=4*n+1 .and. x<4*(m-1)*n+1) then
          h(x:x+4*n-1,x+4*n:x+8*n-1)=r0
          h(x:x+4*n-1,x-4*n:x-1)=l0
        else
          h(x:x+4*n-1,1:4*n)=r0
          h(x:x+4*n-1,x-4*n:x-1)=l0
        end if
  
  end do


   end subroutine hmatr
 end program main

回复列表 (共11个回复)

沙发

Ef(i)=2.6-0.01*i
在循环进行前线分配好Ef(:), 在循环里面你只分配了Ef(i), Ef(i+1)为被初始化

   if(Ef(i)<=eigenvalue(j)<Ef(i+1)) then
     ds(i)=ds(i)+1
   end if
可以考虑改为
   if((Ef(i)<=eigenvalue(j)) .AND. (eigenvalue(j)<Ef(i+1))) then
     ds(i)=ds(i)+1
     exit
   end if

板凳

我再考虑下,根据你说的再试下

3 楼

刚才试了下,运行的时候出错了:
Compiling Fortran...
E:\vfstufdy\taimidu\Heig001.f90
E:\vfstufdy\taimidu\Heig001.f90(75) : Error: Syntax error, found IDENTIFIER 'DS' when expecting one of: ( * :: , <END-OF-STATEMENT> ; : ) + . - % (/ [ ] /) . ' ** / > ...
     write(10,*) Ef(i)  ds(i)
------------------------^
Error executing df.exe.

Heig001.obj - 1 error(s), 0 warning(s)
不知道是怎么回事

4 楼

程序如下:
program main
  use IMSL
!  use linear_operators
  implicit none
  real, parameter::pi=3.14159265
  complex, parameter:: fi=(0.0, 1.0)
  real t,f,w,temp

  integer,dimension(:)::ds(620)
  real,dimension(:)::Ef(620)
  complex*8,allocatable::h(:,:)
  complex*8,allocatable::h0(:,:)
  complex*8,allocatable::r0(:,:)
  complex*8,allocatable::l0(:,:)
  
  complex*8,allocatable::eigenvector(:,:)


  real*8, allocatable:: eigenvalue(:)  !特征值的矩阵
  integer ::i,n,m,j,p
  open(1,file='hmat.txt')
  n=4
  m=4 
  
  allocate(h(4*m*n,4*m*n))
  allocate(h0(4*n,4*n),r0(4*n,4*n),l0(4*n,4*n))
  allocate(eigenvalue(4*n*m)) 
  allocate(eigenvector(4*n*m,4*m*n))

  w=0.0
  t=1.0
  f=0.25

  call hmatr(h, n,m)  
 do  i=1,4*n*m
   do j=1,4*n*m

  write(1,*) h(j,i)

   end do
 end do
close(1)


  eigenvalue=eig(h,w=eigenvector) 
  
  open(10,file='eeig.txt')
  do j=1,4*n*m-1
     p=j
    do i=j+1,4*m*n
       if(eigenvalue(i)<eigenvalue(p)) p=i
    end do
     temp=eigenvalue(j)
     eigenvalue(j)=eigenvalue(p)
     eigenvalue(p)=temp
  end do
  
  do i=1,4*m*n
 
   write(10,*) 'eigenvalue=',eigenvalue(i)
   end do
  close(10)

 
 ds=0
 do j=1,4*m*n
  do i=1,520
    Ef(i)=-2.61+0.01*i
   if((Ef(i)<=eigenvalue(j)) .and. (eigenvalue(j)<Ef(i+1))) then
   ds(i)=ds(i)+1
    exit
   end if
   
   open(10,file='taimidu.txt')
     write(10,*) Ef(i)  ds(i)
   close(10)
   end do
 end do


 stop
 contains
   subroutine hmatr(h,n,m) 
   integer n,m
   complex h(4*m*n,4*m*n)
   complex h0(4*n,4*n)
   complex r0(4*n,4*n)
   complex l0(4*n,4*n)

   integer:: j,n0,x
   complex:: c(4,4),a(4,4),b(4,4),r(4,4),l(4,4)
 
   h0=0.0
   a=0.0
   b=0.0
   h=0.0
   r0=0.0
   l0=0.0
   a(1, 4)=t
   b=transpose(a)
   c=0.0
   r=0.0
   l=0.0

   
      n0=1
   do while(n0<=2*n-1)

      r(2,1)=t*exp(fi*(pi-2.0*pi*f/3.0*(3.0*n0/2.0+0.25)))
      r(3,4)=t*exp(fi*(pi-2.0*pi*f/3.0*(3.0*(n0+1)/2.0+0.25)))
      l(1,2)=t*exp(-fi*(pi-2.0*pi*f/3.0*(3.0*n0/2.0+0.25)))!conjg(r0(2,1))
      l(4,3)=t*exp(-fi*(pi-2.0*pi*f/3.0*(3.0*(n0+1)/2.0+0.25)))!conjg(r0(3,4))


      r0(2*n0,2*n0-1)=r(2,1)
      r0(2*n0+1,2*n0+2)=r(3,4)
      l0(2*n0-1,2*n0)=l(1,2)
      l0(2*n0+2,2*n0+1)=l(4,3)


      c(1,1)=w
      c(2, 1) = t*exp(fi*(-pi+2.0*pi/3.0*f*(1.5*n0+0.25))) 
      c(1, 2) = t*exp(-fi*(-pi+2.0*pi/3.0*f*(1.5*n0+0.25)))
      
      c(2, 2) = w
      c(2, 3) = t
      c(3, 2) = t
      c(3, 3) = w
      c(3, 4) = t*exp(fi*(-pi+2.0*pi/3.0*f*(1.5*(n0+1)+0.25)))
      c(4, 3) = t*exp(-fi*(-pi+2.0*pi/3.0*f*(1.5*(n0+1)+0.25)))
      c(4, 4) = w
    

       h0(2*n0-1,2*n0-1)=c(1,1)
       h0(2*n0,2*n0-1)=c(2,1)
       h0(2*n0-1,2*n0)=c(1,2)
       h0(2*n0,2*n0)=c(2,2)
       h0(2*n0,2*n0+1)=c(2,3)
       h0(2*n0+1,2*n0)=c(3,2)
       h0(2*n0+1,2*n0+1)=c(3,3)
       h0(2*n0+1,2*n0+2)=c(3,4)
       h0(2*n0+2,2*n0+1)=c(4,3)
       h0(2*n0+2,2*n0+2)=c(4,4)   

       if(2*n0-1==1) then

        h0(2*n0-1:2*n0+2,2*n0+3:2*n0+6)=b
        h0(2*n0-1:2*n0+2,4*n-3:4*n)=a

         else if(2*n0-1<4*n-3)  then

          h0(2*n0-1:2*n0+2,2*n0+3:2*n0+6)=b
          h0(2*n0-1:2*n0+2,2*n0-5:2*n0-2)=a

         else 

            h0(4*n-3:4*n,1:4)=b
            h0(2*n0-1:2*n0+2,2*n0-5:2*n0-2)=a
       end if

       n0=n0+2


   end do

 do x=1,4*(m-1)*n+1,4*n
     h(x:x+4*n-1,x:x+4*n-1)=h0
        if( x==1) then
          h(x:x+4*n-1,x+4*n:8*n)=r0
          h(x:x+4*n-1,(m-1)*4*n+1:m*4*n)=l0
        elseif( x>=4*n+1 .and. x<4*(m-1)*n+1) then
          h(x:x+4*n-1,x+4*n:x+8*n-1)=r0
          h(x:x+4*n-1,x-4*n:x-1)=l0
        else
          h(x:x+4*n-1,1:4*n)=r0
          h(x:x+4*n-1,x-4*n:x-1)=l0
        end if
  
  end do


   end subroutine hmatr
 end program main

5 楼

[quote]刚才试了下,运行的时候出错了:
Compiling Fortran...
E:\vfstufdy\taimidu\Heig001.f90
E:\vfstufdy\taimidu\Heig001.f90(75) : Error: Syntax error, found IDENTIFIER 'DS' when expecting one of: ( * :: , <END-OF-STATEMENT> ; : ) + . - % (/ [ ] /) . ' ** / > ...
     write(10,*) Ef(i)  ds(i)
------------------------^
Error executing df.exe.

Heig001.obj - 1 error(s), 0 warning(s)
不知道是怎么回事[/quote]

楼主, Ef(i)  ds(i)之间是不是少了个逗号?? 这个调试能力还是需要具备一些的. 建议还是先弄些简单程序, 锻炼一下.

6 楼

是的,缺少个逗号,但是把逗号加上去还是不能达到预期效果,最后输出数据只有一个,最大的那个数据出来了,从-2.6到2.49这期间的数据都没有,昨天又调试了半天也没什么结果。
ds=0
 do j=1,4*m*n
  do i=1,520
    Ef(i)=-2.61+0.01*i
   if((Ef(i)<=eigenvalue(j)) .and. (eigenvalue(j)<Ef(i+1))) then
   ds(i)=ds(i)+1
    exit
   end if
   
   open(10,file='taimidu.txt')
     write(10,*) Ef(i) , ds(i)
   close(10)
   end do
 end do

7 楼

Ef还没有提取到循环外面. 当你第一次进入j循环的时候Ef(i+1)是未被初始化就被使用了. 试试提出去看看.
ds=0
do i=1,520
  Ef(i)=-2.61+0.01*i
end if
do j=1,4*m*n
  do i=1,520
    if((Ef(i)<=eigenvalue(j)) .and. (eigenvalue(j)<Ef(i+1))) then
      ds(i)=ds(i)+1
      exit
    end if
    open(10,file='taimidu.txt')
    write(10,*) Ef(i) , ds(i)
    close(10)
  end do
end do

8 楼

我再试下。多谢

9 楼

根据你跟我提示的,我修改了程序,也得到了结果,但是有个地方不理想。数据如下:
 -2.602000               1
  -2.602000               2
  -2.602000               3
  -2.602000               4
  -2.594000               1
  -2.594000               2
  -2.594000               3
  -2.594000               4
  -1.874000               1
  -1.874000               2
  -1.874000               3
  -1.874000               4
  -1.810000               1
  -1.810000               2
  -1.810000               3
  -1.810000               4
  -1.418000               1
  -1.418000               2
  -1.418000               3
  -1.418000               4
  -1.306000               1
  -1.306000               2
  -1.306000               3
  -1.306000               4
 -0.3219998               1
 -0.3219998               2
 -0.3219998               3
 -0.3219998               4
 -1.9997712E-03           1
 -1.9997712E-03           2
 -1.9997712E-03           3
 -1.9997712E-03           4
 -1.9997712E-03           5
 -1.9997712E-03           6
 -1.9997712E-03           7
 -1.9997712E-03           8
  0.3100002               1
  0.3100002               2
  0.3100002               3
  0.3100002               4
   1.302000               1
   1.302000               2
   1.302000               3
   1.302000               4
   1.414000               1
   1.414000               2
   1.414000               3
   1.414000               4
   1.806000               1
   1.806000               2
   1.806000               3
   1.806000               4
   1.870000               1
   1.870000               2
   1.870000               3
   1.870000               4
   2.582000               1
   2.582000               2
   2.582000               3
   2.582000               4
   2.590000               1
   2.590000               2
   2.590000               3
   2.590000               4
每个存在数据的地方都给出了1 2 3 4,为什么不直接是4呢?从程序设置上看应该是4啊

ds=0
do i=1,520
  Ef(i)=-2.61+0.01*i
end if
do j=1,4*m*n
  do i=1,520
    if((Ef(i)<=eigenvalue(j)) .and. (eigenvalue(j)<Ef(i+1))) then
      ds(i)=ds(i)+1
      exit
    end if
    open(10,file='taimidu.txt')
    write(10,*) Ef(i) , ds(i)
    
  end do
end do
close(10)  !这个语句调整到这里

10 楼

既然每一个数据只会分配到一个区间里面, 你大可在所有数据都统计完成之后在输出. 

    open(10,file='taimidu.txt')
    write(10,*) Ef(i) , ds(i)
放到最后
open(10,file='taimidu.txt')
do i=1,520
  write(10,*) Ef(i) , ds(i)
end do
close(10)

我来回复

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