时间戳

首页 > 解决方案 > Fortran 中非常少的数字变为负数

fortran - Fortran 中非常少的数字变为负数

问题描述

我正在 Fortran90 中做一个程序，它从i=1给定的i=n地方做一个总和n。总和是sum_{i=1}^{i=n}1/(i*(i+1)*(i+2))。这个和收敛到 0.25。这是代码：

PROGRAM main

INTEGER n(4)
DOUBLE PRECISION s(4)
INTEGER i

  
OPEN(11,FILE='input')
OPEN(12,FILE='output')
DO i=1,4
  READ(11,*) n(i)
END DO

PRINT*,n
CALL suma(n,s)
PRINT*, s

END

 
SUBROUTINE suma(n,s)
INTEGER n(4),j,k
DOUBLE PRECISION s(4),add
s=0 
DO k=1,4
  DO j=1,n(k)
    add=1./(j*(j+1)*(j+2))
    s(k)=s(k)+add 
  END DO 
END DO
  
END SUBROUTINE

输入

178   
1586     
18232 
142705

该output文件现在是空的，我需要对其进行编码。我只是打印结果，它们是： 0.249984481688 0.249999400246 0.248687836759 0.247565846142

问题与变量有关add。当j更大时，add变为负数，并且总和不能很好地收敛。我该如何解决？

标签： fortran

解决方案

这可能是实现您想要的更好的方法。我确实删除了 IO。程序的输出是

% gfortran -o z a.f90 && ./z
   178 0.249984481688392
  1586 0.249999801599584
 18232 0.249999998496064
142705 0.249999999975453

program main

  implicit none  ! Never write a program without this statement

  integer, parameter :: knd = kind(1.d0) ! double precision kind

  integer n(4)
  real(knd) s(4)
  integer i

  n = [178, 1586, 18232, 142705]
  call suma(n, s)
  do i = 1, 4
    print '(I6,F18.15)', n(i), s(i)
  end do

  contains
     !
     ! Recursively, sum a(j+1) = j * a(j) / (j + 1)
     !
     subroutine suma(n, s)
        integer, intent(in) :: n(4)
        real(knd), intent(out) :: s(4)
        real(knd) aj
        integer j, k
        s = 0 
        do k = 1, 4
           aj = 1 / real(1 * 2 * 3, knd)    ! a(1)
           do j = 1, n(k)
              s(k) = s(k) + aj
              aj = j * aj / (j + 3)
           end do 
        end do
     end subroutine

end program main

推荐阅读