Fortran自定义arctan子程序无法正常运行问题咨询

Hey there! Let's work through your Fortran arctan subroutine issue step by step. I see you're trying to implement the Taylor series for arctan to compute π using Machin's formula, and your current subroutine isn't working as expected. Let's break this down and fix it.

First, let's spot the issues in your code snippet #

Your current approach has a few key problems:

• Incorrect loop step: The Taylor series for arctan uses terms with exponents increasing by 2 (1, 3, 5, 7...), but you're using a step of 4, which skips half the necessary terms.
• Incomplete/incorrect term calculation: Your code cuts off at x**(n..., and the logic for grouping terms is flawed.
• Inefficient power calculation: Directly computing x**num for large num wastes resources and can introduce precision loss. A recursive term calculation is better here.

Fixed Arctan Subroutine (Using Taylor Series with Recursion) #

The Taylor series for arctan(x) is:

arctan(x) = x - x³/3 + x⁵/5 - x⁷/7 + ... = Σ (-1)^k * x^(2k+1)/(2k+1) for k ≥ 0

Instead of calculating each power directly, we can compute each term recursively from the previous one to boost efficiency and precision. Here's the corrected subroutine:

subroutine arctan(x, n, arc)
    implicit none
    real*8, intent(in)  :: x       ! Input value for arctan
    integer, intent(in)  :: n       ! Number of terms to compute (more = higher precision)
    real*8, intent(out) :: arc      ! Output arctan(x) value
    integer :: k
    real*8 :: term, x_squared

    ! Initialize variables
    arc = 0.0d0
    x_squared = x * x
    term = x  ! First term (k=0: x^1 / 1)

    ! Compute each term recursively
    do k = 1, n
        arc = arc + term
        ! Update term for next iteration: sign flips, exponent increases by 2, denominator increases by 2
        term = term * (-x_squared) * (2*k - 1) / (2*k + 1)
    end do

end subroutine arctan

Full Program to Compute π with Machin's Formula #

Here's a complete main program that uses the fixed subroutine to calculate π via π = 16*arctan(1/5) - 4*arctan(1/239):

program compute_pi
    implicit none
    real*8 :: pi, atan_1_5, atan_1_239
    integer :: num_terms

    ! 10 terms are enough for ~15 decimal digits of precision
    num_terms = 10

    ! Calculate the two arctan values
    call arctan(1.0d0/5.0d0, num_terms, atan_1_5)
    call arctan(1.0d0/239.0d0, num_terms, atan_1_239)

    ! Apply Machin's formula
    pi = 16.0d0 * atan_1_5 - 4.0d0 * atan_1_239

    ! Print results
    print *, "Computed π with ", num_terms, " terms:"
    print *, pi
    print *, "Difference from actual π: ", pi - 4.0d0*datan(1.0d0)

end program compute_pi

Key Notes #

• implicit none: Always include this to avoid accidental implicit type declarations (a common Fortran pitfall for beginners).
• Recursive term calculation: By updating each term from the previous one, we avoid expensive high-power computations and keep the numerical calculation stable.
• Precision: Using real*8 (double-precision) is critical for getting accurate π values—single-precision (real*4) won't be sufficient for this use case.

内容的提问来源于stack exchange，提问作者Arashi Nakamura