因此,作为一项任务,我的任务是编写一个函数,当给定一个 x 时,它会从中计算出相应的一阶贝塞尔函数。公式如下:https ://youtu.be/vBOYr3m2M8E?t=48(抱歉没有足够的声誉来发布照片)。尽管我的条件是,当第 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