时间戳

首页 > 解决方案 > Fortran arctan 子例程未按预期工作

fortran - Fortran arctan 子例程未按预期工作

问题描述

我通常是 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 都会重置，所以我知道这不是问题。我认为问题在于我在声明firstand之前如何调用子例程second，但我不知道我能做些什么来改进它们。

我究竟做错了什么？

编辑：我实际上解决了这个问题，而实际问题只是 fortran 决定它不想这样做approx = first - second 并且正在制作它，所以大约 == second 我不知道为什么，但我通过用以下语句替换该语句解决了这个问题：

approx = (second-first)
approx = approx *(-1)

尽管看起来很愚蠢，但它现在可以完美运行，结果为 3.1415926535897940！

标签： fortransubroutine

解决方案

该问题是由 arctan 调用中变量 arc 的不同类型（单精度/双精度）和子程序的实现引起的。10000 次的迭代次数...太多了，可能会导致数值问题，仅 100 次就足够了（而且速度更快...）。

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

只需在主程序中将其设为双精度，您就会得到所需的答案。

推荐阅读