【问题标题】:How do I make all the calculations in double precision in fortran?如何在 fortran 中以双精度进行所有计算?
【发布时间】:2015-04-27 19:04:11
【问题描述】:

在下面给出的 Fortran 代码中,我将所有涉及计算 PI 的数字都设为双精度,但我得到的 PI 值只是一个实数,末尾有大量的零或 9。如何让程序以双精度给出 PI?我正在使用gfortran 编译器。

  !This program determines the value of pi using Monte-Carlo algorithm.
  program findpi
  implicit none
  double precision :: x,y,radius,truepi,cnt
  double precision,allocatable,dimension(:) :: pi,errpi
  integer :: seedsize,i,t,iter,j,k,n
  integer,allocatable,dimension(:) :: seed

  !Determining the true value of pi to compare with the calculated value
  truepi=4.D0*ATAN(1.D0)

  call random_seed(size=seedsize)
  allocate(seed(seedsize))
  do i=1,seedsize
     call system_clock(t) !Using system clock to randomise the seed to 
                          !random number generator
     seed(i)=t
  enddo
  call random_seed(put=seed)

  n=2000         !Number of times value of pi is determined
  allocate(pi(n),errpi(n))
  do j=1,n
     iter=n*100  !Number of random points
     cnt=0.D0
     do i=1,iter
        call random_number(x)
        call random_number(y)
        radius=sqrt(x*x + y*y)
        if (radius < 1) then
           cnt = cnt+1.D0
        endif
     enddo
     pi(j)=(4.D0*cnt)/dble(iter)
     print*, j,pi(j)
  enddo

  open(10,file="pi.dat",status="replace")
  write(10,"(F15.8,I10)") (pi(k),k,k=1,n)

  call system("gnuplot --persist piplot.gnuplot")

end program findpi

【问题讨论】:

    标签: fortran fortran90 gfortran


    【解决方案1】:

    您的计算是双精度的,但我发现两个问题:

    • 第一个是系统错误...你通过
    • 确定pi
    pi(j)=(4.D0*cnt)/dble(iter)
    

    iter 最多为 2000*100,因此1/iter 至少为 5e-6,因此您无法解析任何查找器;-)

    • 第二个问题是您的 IO 例程以单精度打印结果!行
    write(10,"(F15.8,I10)") (pi(k),k,k=1,n)
    

    更具体地说,需要调整格式说明符"(F15.8,I10)"。目前它告诉编译器总共使用 15 个字符来打印数字,小数点后有 8 位。作为第一个措施,您可以使用*:

    write(10,*) (pi(k),k,k=1,n)
    

    这总共使用 22 个字符,所有 15 位数字用于双精度:

    write(10,"(F22.15,I10)") (pi(k),k,k=1,n)
    

    【讨论】:

      猜你喜欢
      • 2015-06-06
      • 2010-12-04
      • 1970-01-01
      • 2017-05-18
      • 1970-01-01
      • 2017-01-06
      • 1970-01-01
      • 1970-01-01
      • 1970-01-01
      相关资源
      最近更新 更多