UW AMath High Performance Scientific Computing
 
AMath 483/583 Class Notes
 
Spring Quarter, 2013

Previous topic

Fortran subroutines and functions

Next topic

Array storage in Fortran

This Page

Fortran examples: Taylor seriesΒΆ

Here is an example code that uses the Taylor series for \(exp(x)\) to estimate the value of this function at \(x=1\):

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
! $UWHPSC/codes/fortran/taylor.f90

program taylor

    implicit none                  
    real (kind=8) :: x, exp_true, y
    real (kind=8), external :: exptaylor
    integer :: n

    n = 20               ! number of terms to use
    x = 1.0
    exp_true = exp(x)
    y = exptaylor(x,n)   ! uses function below
    print *, "x = ",x
    print *, "exp_true  = ",exp_true
    print *, "exptaylor = ",y
    print *, "error     = ",y - exp_true

end program taylor

!==========================
function exptaylor(x,n)
!==========================
    implicit none

    ! function arguments:
    real (kind=8), intent(in) :: x
    integer, intent(in) :: n
    real (kind=8) :: exptaylor

    ! local variables:
    real (kind=8) :: term, partial_sum
    integer :: j

    term = 1.
    partial_sum = term

    do j=1,n
        ! j'th term is  x**j / j!  which is the previous term times x/j:
        term = term*x/j   
        ! add this term to the partial sum:
        partial_sum = partial_sum + term   
        enddo
     exptaylor = partial_sum  ! this is the value returned
end function exptaylor

Running this code gives:

x =    1.00000000000000
exp_true  =    2.71828182845905
exptaylor =    2.71828174591064
error     =  -8.254840055954560E-008

Here’s a modification where the number of terms to use is computed based on the size of the next term in the series. The function has been rewritten as a subroutine so the number of terms can be returned as well.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
! $UWHPSC/codes/fortran/taylor_converge.f90

program taylor_converge

    implicit none                  
    real (kind=8) :: x, exp_true, y, relative_error
    integer :: nmax, nterms, j

    nmax = 100

    print *, "     x         true              approximate          error         nterms"
    do j = -20,20,4
       x = float(j)                      ! convert to a real
       call exptaylor(x,nmax,y,nterms)   ! defined below
       exp_true = exp(x)
       relative_error = abs(y-exp_true) / exp_true
       print '(f10.3,3d19.10,i6)', x, exp_true, y, relative_error, nterms
       enddo

end program taylor_converge

!====================================
subroutine exptaylor(x,nmax,y,nterms)
!====================================
    implicit none

    ! subroutine arguments:
    real (kind=8), intent(in) :: x
    integer, intent(in) :: nmax
    real (kind=8), intent(out) :: y
    integer, intent(out) :: nterms

    ! local variables:
    real (kind=8) :: term, partial_sum
    integer :: j

    term = 1.
    partial_sum = term

    do j=1,nmax
        ! j'th term is  x**j / j!  which is the previous term times x/j:
        term = term*x/j   
        ! add this term to the partial sum:
        partial_sum = partial_sum + term   
        if (abs(term) < 1.d-16*partial_sum) exit
        enddo
     nterms = j       ! number of terms used
     y = partial_sum  ! this is the value returned
end subroutine exptaylor

This produces:

   x         true              approximate          error         nterms
-20.000   0.2061153622D-08   0.5621884472D-08   0.1727542678D+01    95
-16.000   0.1125351747D-06   0.1125418051D-06   0.5891819580D-04    81
-12.000   0.6144212353D-05   0.6144213318D-05   0.1569943213D-06    66
 -8.000   0.3354626279D-03   0.3354626279D-03   0.6586251980D-11    51
 -4.000   0.1831563889D-01   0.1831563889D-01   0.1723771005D-13    34
  0.000   0.1000000000D+01   0.1000000000D+01   0.0000000000D+00     1
  4.000   0.5459815003D+02   0.5459815003D+02   0.5205617665D-15    30
  8.000   0.2980957987D+04   0.2980957987D+04   0.1525507414D-15    42
 12.000   0.1627547914D+06   0.1627547914D+06   0.3576402292D-15    51
 16.000   0.8886110521D+07   0.8886110521D+07   0.0000000000D+00    59
 20.000   0.4851651954D+09   0.4851651954D+09   0.1228543295D-15    67

Comments:

  • Note the use of exit to break out of the loop.
  • Note that it is getting full machine precision for positive values of x but, as expected, the accuracy suffers for negative x due to cancellation.