[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
|