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

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


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

fcode 发表于 2016-1-21 17:30:42

http://v.fcode.cn/video-debugger.html
请查看此视频教程,学习如何debug调试

pasuka 发表于 2016-1-21 18:47:47

fcode 发表于 2016-1-21 17:30
http://v.fcode.cn/video-debugger.html
请查看此视频教程,学习如何debug调试

马上开始放寒假了,这类帖子也多起来了。。。
页: [1]
查看完整版本: 求助各位大神帮忙看看程序哪里有问题,解出来X为无穷大...