!___________________________________________________________________________________________________ ! 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 !