求助各位大神帮忙看看程序哪里有问题，解出来X为无穷大...

sophierunrun · 发表于 2016-1-21 15:29:04

[Fortran] 纯文本查看 复制代码

module typedef
   type distribution
      real(kind=8) :: data1
      real(kind=8) :: data2
   end type
end module

module numerical
  implicit none
  real*8,parameter :: eps=0.000001
contains
  real function newton(a,f,df)
  implicit none
  real(kind=8) :: a
  real, external :: f
  real, external :: df
  real(kind=8) :: b, fb
  integer :: t=0

  b=a-f(a)/df(a)
  fb=f(b)
  do while(abs(fb)>eps)
    a=b
    b=a-f(a)/df(a)
    fb=f(b)
        t=t+1
  end do
  newton=b
  return
  end function newton

  real function func(x)
  implicit none
  real, parameter :: n=50.0
  real(kind=8) :: x
  real(kind=8):: n1,n2
  common n1,n2
  func=n*n1*log(x)-n*n2-n*x
  return
  end function func

  real function dfunc(x)
  implicit none
  real, parameter :: n=50.0
  real(kind=8):: o1,o2
  common o1,o2
  real(kind=8) :: x
  dfunc=n*o1/x-n
  return
  end function dfunc
end module

program main
  use typedef
  use numerical
  implicit none
  integer, parameter :: n=50
  integer, parameter :: fileid=10
  type(distribution) :: p(n)
  type(distribution) :: total
  integer :: i, error
  real(kind=8):: m1,m2
  common m1,m2
  real(kind=8) :: x0
  real(kind=8) :: ans

  open(fileid,file="data.txt",status="old",iostat=error)
        if(error/=0) then
          write(*,*) "open data.txt fail"
        stop
        end if
  total=distribution(0,0)
  do i=1,n
         read(fileid,*) p(i)%data1, p(i)%data2
         total%data1=total%data1+p(i)%data1
         total%data2=total%data2+p(i)%data2
  end do
   m1=total%data1/real(n)
   m2=total%data2/real(n) 

   x0=m1
   ans=newton(x0,func,dfunc)
   write(*,"('x=',f8.4)") ans
stop
end program