0

So as an assignment I was given the task to write a function that when given an x, calculates the corresponding first order Bessel function from it. The equation is as follows: https://youtu.be/vBOYr3m2M8E?t=48 (sorry don't have enough reputation to post a photo). My implementation goes on infinitely despite the fact my condition, which is when the r-th summation value is less than some epsilon (the do-while code), mathematically should eventually fail (because as n approaches infinity, n!(n+1)! >> (x/2)^n). I've traced out the input that I can by pausing after execution and I noticed after about the 5th iteration that my program calculates an incorrect value (-67 instead of 40) but I'm confused why this happens, especially since it works initially. I've also searched online for examples so I am aware of the presence of a library that does this for me, but that defeats the purpose of the assignment. I was hoping someone could point out why this is occurring and maybe also let me know if my implementation is incorrect in any other aspects.

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
Ian Bush
  • 6,996
  • 1
  • 21
  • 27
badejayo
  • 15
  • 5
  • Too late to look properly now but those expressions involving factorials are probably overflowing the integers you are using for them. Also note a) Real( 8 ) is not portable and so poor practice b) bessel_j1 is part of standard Fortran, which might be useful to test your code c) Please do NOT post pictures here, write the formula as best you can – Ian Bush Mar 28 '21 at 21:00
  • Bessel functions are intrinsically part of Fortran language. You do not need to implement them. I understand this is a HW though. Here are the intrinsic function names: https://gcc.gnu.org/onlinedocs/gcc-4.7.4/gfortran/BESSEL_005fJN.html https://gcc.gnu.org/onlinedocs/gfortran/BESSEL_005fYN.html – Scientist Mar 28 '21 at 21:02
  • Also, you do not need to implement factorial. It exists in Fortran as `log_gamma(x)` which returns the natural log of the factorial of (x-1). This also avoids the overflow mentioned by @IanBush. To get the factorial instead of the log, try `gamma()`: https://gcc.gnu.org/onlinedocs/gcc-4.9.4/gfortran/GAMMA.html#GAMMA – Scientist Mar 28 '21 at 21:05
  • 1
    Which integer order Bessel function are you trying to compute? Actually, doesn't matter as it appears you're trying to use the infinite series solution. For x > 2, you're going to suffer from catastrophic cancellation. You also do not want to compute a factorial. You can write the (n+1)th term in terms of the nth term. With that simple relation, you can recursively evaluate the (n+1)th. – steve Mar 28 '21 at 21:18
  • I'd like to thank you all for taking the time to answer the question. I found the gamma solutions to be very helpful, although they didn't solve the problem of the infinite looping. I also will take into account your hint for the Bessel function in regards to testing. – badejayo Mar 28 '21 at 22:30
  • Steve, in regards to your comment, I was primarily testing x = 10 in hopes of getting 0.04 but I seemed to continually get infinite looping. When you said I could write the (n+1)th term in terms of the nth term was that for calculating factorials or calculating the bessel function? I'm just asking because my teacher gave a similar hint for the bessel equation so I wondered if you meant the same thing. – badejayo Mar 28 '21 at 22:33
  • So, you're compute an infinite sum: J0(x) = a0 + a1 + a2 + .... Let z = x / 2. The 1st term is a0 = 1. The 2nd term is a1 = - z * z. By induction, one can show that a(n+1) = (- z * z / (n * n)) * an for n = 1, 2, 3 .... You can recursively compute the an. This is, however, a poor way to compute J0(x). Simply inspect the magnitudes of the an shows the cancellation problem. – steve Mar 28 '21 at 23:48

1 Answers1

0

Here is a proof-of-concept for summing the infinite series for J0(x) in default real. You are welcomed to use for your assignment if you can explain the code; especially, lines I did not comment.

Caution: Do not compile this code with gfortran's -ffast-math option.

!
! 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
steve
  • 657
  • 1
  • 4
  • 7