24小时热门版块排行榜

返回列表

当前主题已经存档。

lvjian8596

木虫 (小有名气)

应助: 0 (幼儿园)
贵宾: 0.014
金币: 3531.7
红花: 2
帖子: 296
在线: 209.1小时
虫号: 545009
注册: 2008-04-13
性别: GG
专业: 凝聚态物性I:结构、力学和

[交流] 【求助】关于在对数网格上进行SIMPSON积分的问题

我想用fortran程序编写一个关于积分的程序
首先需要把函数产生在一个对数网格上，然后在这个对数网格上进行simpson积分！
看了PWSCF上的一段关于产生对数网格的程序，不是很理解，希望大家帮忙！
program main
  implicit none
  integer    :: i,j,k
  real(kind=8) :: xmin,zmesh,dx,mesh,rab(1000),rab1(1000),r(1000),x

      mesh  = 1000
      xmin  = -5.0d0
      zmesh = 6.0d0
      dx =  0.025d0
      DO i = 1, mesh
         x = xmin + DBLE(i-1) * dx
         r(i)  = exp(x) / zmesh
         IF( r(i) > 1000.0d0 ) EXIT
      END DO
      rab  = dx * r
end
simpson 积分的子程序是：
subroutine simpson(mesh,func,rab,asum)
  !-----------------------------------------------------------------------
  !    simpson's rule integrator for function stored on the
  !    radial logarithmic mesh

  implicit none

  integer :: i, mesh
  real(8) ::  rab(mesh), func(mesh), f1, f2, f3, r12, asum
  asum = 0.0d0
  r12 = 1.0d0 / 12.0d0
  f3  = func(1) * rab(1) * r12

  do i = 2,mesh-1,2
   f1 = f3
   f2 = func(i) * rab(i) * r12
   f3 = func(i+1) * rab(i+1) * r12
   asum = asum + 4.0d0*f1 + 16.0d0*f2 + 4.0d0*f3
  enddo
  return
end subroutine simpson
我不理解的地方是：
r(i)应该是网格上的点
rab(i)应该是产生的网格之间的距离吧???
那为什么这个距离是 rab  = dx * r
而不是 rab(i)=r(i+1)-r(i)

请大家帮助！

[ Last edited by woshilsh on 2008-12-14 at 20:40 ]

回复此楼