【问题标题】:Fortran arctan subroutine is not working as expectedFortran arctan 子例程未按预期工作
【发布时间】: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 都会重置,所以我知道这不是问题。我认为问题在于我在声明 firstsecond 之前如何调用子例程,但我不知道我可以做些什么来改进它们。

我做错了什么?

编辑: 我实际上解决了问题,实际问题只是 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


【解决方案1】:

问题是由不同类型(单/双精度)的 变量arc在arctan的调用和子程序的执行中。 10000 次的迭代次数...太多了,可能会导致数值问题,仅 100 次就足够了(而且速度更快...)。

提示:始终对所有 progs 和过程使用隐式 none。在这里,编译器会立即告诉您您忘记声明 arc...

只要在主程序中设置双精度,你就会得到想要的答案。

【讨论】:

  • 感谢您的意见。我确保将 arc 更改为双精度,并在我的程序和子例程中添加了隐式 none。我还将我的 n 更改为 100。不幸的是,我仍然得到 1.67363040082988898E-002 作为我的答案。
  • @ArashiNakamura 请提供显示问题的最小完整代码,即带有被调用子程序并且我们可以运行的小主程序。
猜你喜欢
  • 1970-01-01
  • 2018-04-27
  • 2018-03-21
  • 1970-01-01
  • 1970-01-01
  • 2022-01-13
  • 1970-01-01
  • 1970-01-01
  • 2012-07-29
相关资源
最近更新 更多