!___________________________________________________________________________________________________
! Aleksandar Donev
! PHY201 Solutions to Worksheet
! December, 2000, MSU
!___________________________________________________________________________________________________
PROGRAM Erf_Series
!
IMPLICIT NONE
INTEGER, PARAMETER :: sp = KIND (0.0E0), dp = KIND (0.0D0)
INTEGER, PARAMETER :: wp = dp ! Use single precision in this run
INTEGER, PARAMETER :: max_iterations = 100
INTEGER :: status ! I/O status
INTEGER :: k ! Iteration counter
REAL (KIND=wp) :: x, erf ! x and erf(x)
REAL (KIND=wp) :: sum, factor, expression ! A few auxiliary variables
REAL (KIND=wp) :: error, epsilon ! Error allowed in erf(x) and in the sum
REAL (KIND=wp) :: pi ! The constant pi
!
MainLoop: DO
!
WRITE (UNIT=*, FMT="(A)", ADVANCE="NO") "Enter the value of x and epsilon ('S' to stop): "
READ (UNIT=*, FMT=*, IOSTAT=STATUS) x, error
IF (status /= 0) EXIT MainLoop
!
pi = 4 * Atan (1.0_wp)! This is a good way of calculating pi to the needed precision
factor = 2 * x / Sqrt (pi)! The factor in front of the sum
epsilon = Abs (error/factor)! This is the error allowed *in* the sum
!
Magnitude: IF (Abs(x) < 5.0_wp) THEN ! We can use the series provided
sum = 1.0_wp ! This is E(0)
expression = 1.0_wp ! This is E(0) as well
k = 0 ! Initialize counter
Summation: DO
k = k + 1 ! Iteration counter
! Calculate the next term to be added
expression = - REAL (2*k-1, wp) / REAL (2*k+1, wp) / REAL (k, wp) * (x**2) * expression
! Now check for convergence:
IF (Abs(expression) < epsilon) THEN ! Converged
WRITE (UNIT=*, FMT=*) "Convergence was achieved with ", k, " terms in the sum"
erf = factor * sum
WRITE (UNIT=*, FMT=*) "Erf (", x, ") = ", erf
EXIT Summation
ELSE IF (k >= max_iterations) THEN ! Not converged
WRITE (UNIT=*, FMT=*) "Convergence was not achieved!"
EXIT Summation
ELSE ! Keep on adding terms
sum = sum + expression ! Add the term to the sum
CYCLE Summation ! This can be ommitted
END IF
END DO Summation
ELSE ! If x is large in magnitude, we can not use the above series easily
erf = 1.0_wp
WRITE (UNIT=*, FMT=*) "The value of x is too large! An approximate answer is:"
WRITE (UNIT=*, FMT=*) "Erf (", x, ") = ", erf
END IF Magnitude
END DO MainLoop
!
END PROGRAM Erf_Series
!