|
High Performance Scientific Computing
AMath 483/583 Class Notes Spring Quarter, 2013 |
|
|
High Performance Scientific Computing
AMath 483/583 Class Notes Spring Quarter, 2013 |
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.