【发布时间】:2018-04-29 06:30:04
【问题描述】:
我通常是 Fortran 新手,我有一个项目,我的教授希望全班尝试找到 pi。为此,他希望我们创建自己的 arctan 子程序并使用这个特定的公式:pi = 16*arctan(1/5) - 4*arctan(1/239)。
因为教授不让我使用内置的ATAN函数,所以我做了一个近似它的子程序:
subroutine arctan(x,n,arc)
real*8::x, arc
integer::n, i
real*8::num, nm2
arc = 0.0
do i=1,n,4
num = i
nm2 = num+2
arc = arc+((x**num)/(num)) - (x**(nm2)/(nm2))
enddo
end subroutine arctan
此子例程基于用于反正切近似的泰勒级数,并且似乎可以完美运行,因为我通过调用它对其进行了测试。
real*8:: arc=0.0, approx
call arctan(1.d0,10000000,arc)
approx = arc*4
我从我的主程序中调用了这个,它应该返回 pi,我得到了
approx = 3.1415926335902506
这对我来说已经足够接近了。当我尝试做时出现问题
pi = 16*arctan(1/5) - 4*arctan(1/239)。我试过这个:
real*8:: first, second
integer:: n=100
call arctan((1.d0/5.d0), n, arc)
first = 16*arc
call arctan((1.d0/239.d0), n, arc)
second = 4.d0*arc
approx = first - second
不知何故approx = 1.67363040082988898E-002,这显然不是pi。每次调用 arctan 子例程时 arc 都会重置,所以我知道这不是问题。我认为问题在于我在声明 first 和 second 之前如何调用子例程,但我不知道我可以做些什么来改进它们。
我做错了什么?
编辑:
我实际上解决了问题,实际问题只是 fortran 决定它不想做approx = first - second
并且正在使它成为大约 == 秒我不知道为什么,但我通过用以下语句替换该语句解决了这个问题:
approx = (second-first)
approx = approx *(-1)
尽管看起来很愚蠢,但它现在可以完美运行,结果为 3.1415926535897940!
【问题讨论】:
-
您的问题可能与您将迭代次数
n传递给 arctan 函数的方式有关。首先,调用者不应该猜测泰勒级数中需要多少项才能收敛。但是,其次,您可能会发现,将十亿 (!) 项相加,实际上引入的舍入误差要比在数量不多的项之后终止该系列要多得多。如果您只需要计算 arctan(1/5)(或更小),那么该级数应该会很快收敛。 -
在“错误的示例中,您声明了 arc?使其成为可以编译和运行的 MWE(最小工作示例)。
-
我认为在这种情况下,代码中的
IMPLICIT NONE或命令行中的-fimplicit-none可能是你的朋友。 -
如果您正在学习 Fortran,请不要使用 real*8 - 它不是 Fortran 的一部分,也从未成为 Fortran 的一部分。而是了解 Fortran 种类机制,这是应该如何完成的。
-
我非常怀疑您的“修复”,这将是一个主要的编译器错误。由于您没有发布您的代码,我们无法看到真正出了什么问题(以及您可能已经更改了更多,而您没有记住它/没有写在您的编辑中)。
标签: fortran subroutine