Fortran Coder

查看: 8652|回复: 2
打印 上一主题 下一主题

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

[复制链接]

2

帖子

2

主题

0

精华

新人

F 币
17 元
贡献
9 点
跳转到指定楼层
楼主
发表于 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


分享到:  微信微信
收藏收藏 点赞点赞 点踩点踩

2015

帖子

12

主题

5

精华

论坛跑堂

臭石头雪球

F 币
1574 元
贡献
676 点

美女勋章热心勋章星光勋章新人勋章贡献勋章管理勋章帅哥勋章爱心勋章规矩勋章元老勋章水王勋章

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

490

帖子

4

主题

0

精华

大宗师

F 币
3298 元
贡献
1948 点

水王勋章元老勋章热心勋章

板凳
发表于 2016-1-21 18:47:47 | 只看该作者
fcode 发表于 2016-1-21 17:30
http://v.fcode.cn/video-debugger.html
请查看此视频教程,学习如何debug调试

马上开始放寒假了,这类帖子也多起来了。。。
您需要登录后才可以回帖 登录 | 极速注册

本版积分规则

捐赠本站|Archiver|关于我们 About Us|小黑屋|Fcode ( 京ICP备18005632-2号 )

GMT+8, 2024-11-13 14:24

Powered by Tencent X3.4

© 2013-2024 Tencent

快速回复 返回顶部 返回列表