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
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