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

#### Previous topic

Fortran subroutines and functions

#### Next topic

Array storage in Fortran

# 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