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