在 Fortran 函数中重新实现 Bessel 函数导致无限循环

Question

因此，作为一项任务，我的任务是编写一个函数，当给定一个 x 时，它会从中计算相应的一阶贝塞尔函数。等式如下：https://youtu.be/vBOYr3m2M8E?t=48（抱歉没有足够的声誉来 post 一张照片）。尽管我的条件是，当第 r 个求和值小于某个 epsilon（do-while 代码）时，我的实现仍在无限进行，数学上最终应该失败（因为当 n 接近无穷大时，n!(n+1 )!>> (x/2)^n).我已经通过在执行后暂停来追踪我可以输入的内容，并且在大约第 5 次迭代后我注意到我的程序计算了一个不正确的值（-67 而不是 40）但是我很困惑为什么会发生这种情况，特别是因为它最初有效.我还在网上搜索了示例，所以我知道有一个图书馆可以为我做这件事，但这违背了作业的目的。我希望有人能指出为什么会发生这种情况，如果我的实施在任何其他方面不正确，也可能会告诉我。

implicit none

real (kind = 8) :: x, eps, current, numerator, iteration
integer :: counter, m, denominator

eps = 1.E-3  
counter = 0
m = 1

print*, 'What is your x value? '  
read*, x

current = 1/factorial(m)
print*, current

if (abs(((x / 2) ** m) * current) < eps) THEN
    counter = 1 
    current = ((x / 2) ** m) * current
    print*, current
else 
counter = 1
iteration = current
do while(abs(iteration) > eps)
    numerator = ((-1) ** counter) * ((x / 2) ** (counter * 2))
    denominator = (factorial(counter) * factorial(counter + m))
    iteration = (numerator / denominator)
    current = current + iteration
    counter = counter + 1 
    print*, counter
    print*, current
end do
current = ((x / 2) ** m) * current 
end if 

CONTAINS

recursive function factorial(n) result(f)
integer :: f, n 

if (n == 1 .or. n == 0) THEN
    f = 1 
else 
    f = n * factorial(n - 1)
end if
end function factorial

end program bessel

Answer 1

这里是对 J0(x) 的无穷级数求和的概念验证（默认实数）。如果你能解释代码，欢迎你用于你的任务；特别是，我没有评论的台词。

警告：不要使用 gfortran 的 -ffast-math 选项编译此代码。

!
! Define a named constant for default real.
!
module mytypes

  implicit none  ! Always include this line

  private        ! Make everything private
  public sp      ! Expose only sp

  integer, parameter :: sp = kind(1.e0)  ! Default real

end module mytypes
!
! Compute the zeroth order Bessel of real argument via summation. This
! uses direct summation of the J0(x) = a0 + a1 + a2 + ..., which is not
! a good idea for |x| > 2 due to catastrophic cancellation.  This also
! has really bad results near zeros of J0(x).
! 
module bessel

  use mytypes, only : sp

  implicit none  ! Always include this line

  private        ! Make everything private
  public j0f     ! Expose only j0f

  contains 

     impure elemental function j0f(x) result(res)

        real(sp) res
        real(sp), intent(in) :: x

        integer m, n
        real(sp), volatile, save :: tiny = 1.e-30_sp
        real(sp) a0, c, t, y, z2

        z2 = abs(x)

        if (z2 < scale(1._sp, -digits(z2)/2)) then
           res = 1 - tiny
           return
        end if

        z2 = (z2 / 2)**2
        a0 = 1
        if (x < 2) then
           c = 0
           res = a0
           do m = 1, 5
              a0 = - z2 * a0 / m**2
              y = a0 - c
              t = res + y
              c = (t - res) - y
              res = t
           end do
           if (x < 1) return
           a0 = - z2 * a0 / m**2
           y = a0 - c
           t = res + y
           c = (t - res) - y
           res = t
           m = m + 1
           a0 = - z2 * a0 / m**2
           y = a0 - c
           t = res + y
           c = (t - res) - y
           res = t
        else
           n = 4 * x
           c = 0
           res = a0
           do m = 1, n
              a0 = - z2 * a0 / m**2
              y = a0 - c
              t = res + y
              c = (t - res) - y
              res = t
           end do
        end if
     end function j0f
end module bessel

program foo

  use bessel, only : j0f
  use mytypes, only : sp

  integer, parameter :: n = 100
  integer i

  real(sp) e(n), j(n), x(n)
  real(sp), parameter :: xmax = 10

  x = [(i, i = 0, n - 1)] * (xmax / (n - 1))
  e = bessel_j0(x)
  j = j0f(x)

  do i = 1, n
     write(*,'(3F12.7,ES12.4)') x(i), e(i), j(i), &
     & abs((e(i) - j(i)) / e(i)) / epsilon(1._sp)
  end do

end program foo

在 Fortran 函数中重新实现 Bessel 函数导致无限循环

Reimplementing Bessel Function in Fortran Function Causing Infinite Looping

fortran

bessel-functions