! D05BEF Example Program Text ! Mark 23 Release. NAG Copyright 2011. MODULE d05befe_mod ! D05BEF Example Program Module: ! Parameters and User-defined Routines ! .. Use Statements .. USE nag_library, ONLY : nag_wp ! .. Implicit None Statement .. IMPLICIT NONE ! .. Parameters .. INTEGER, PARAMETER :: iorder = 4 INTEGER, PARAMETER :: nmesh = 2**6 + 2*iorder - 1 INTEGER, PARAMETER :: nout = 6 INTEGER, PARAMETER :: lct = nmesh/32 + 1 INTEGER, PARAMETER :: & lwk = (2*iorder+6)*nmesh + 8*iorder*iorder - 16*iorder + 1 CONTAINS FUNCTION sol1(t) ! .. Use Statements .. USE nag_library, ONLY : x01aaf ! .. Implicit None Statement .. IMPLICIT NONE ! .. Function Return Value .. REAL (KIND=nag_wp) :: sol1 ! .. Scalar Arguments .. REAL (KIND=nag_wp), INTENT (IN) :: t ! .. Local Scalars .. REAL (KIND=nag_wp) :: c, pi, t1 ! .. Intrinsic Functions .. INTRINSIC exp, sqrt ! .. Executable Statements .. t1 = 1.0_nag_wp + t c = 1.0_nag_wp/sqrt(2.0_nag_wp*x01aaf(pi)) sol1 = c*(1.0_nag_wp/(t**1.5_nag_wp))*exp(-t1*t1/(2.0_nag_wp*t)) RETURN END FUNCTION sol1 FUNCTION sol2(t) ! .. Implicit None Statement .. IMPLICIT NONE ! .. Function Return Value .. REAL (KIND=nag_wp) :: sol2 ! .. Scalar Arguments .. REAL (KIND=nag_wp), INTENT (IN) :: t ! .. Executable Statements .. sol2 = 1.0_nag_wp/(1.0_nag_wp+t) RETURN END FUNCTION sol2 FUNCTION ck1(t) ! .. Implicit None Statement .. IMPLICIT NONE ! .. Function Return Value .. REAL (KIND=nag_wp) :: ck1 ! .. Scalar Arguments .. REAL (KIND=nag_wp), INTENT (IN) :: t ! .. Intrinsic Functions .. INTRINSIC exp ! .. Executable Statements .. ck1 = exp(-0.5_nag_wp*t) RETURN END FUNCTION ck1 FUNCTION cf1(t) ! .. Use Statements .. USE nag_library, ONLY : x01aaf ! .. Implicit None Statement .. IMPLICIT NONE ! .. Function Return Value .. REAL (KIND=nag_wp) :: cf1 ! .. Scalar Arguments .. REAL (KIND=nag_wp), INTENT (IN) :: t ! .. Local Scalars .. REAL (KIND=nag_wp) :: a, pi, t1 ! .. Intrinsic Functions .. INTRINSIC exp, sqrt ! .. Executable Statements .. t1 = 1.0_nag_wp + t a = 1.0_nag_wp/sqrt(x01aaf(pi)*t) cf1 = -a*exp(-0.5_nag_wp*t1*t1/t) RETURN END FUNCTION cf1 FUNCTION cg1(s,y) ! .. Implicit None Statement .. IMPLICIT NONE ! .. Function Return Value .. REAL (KIND=nag_wp) :: cg1 ! .. Scalar Arguments .. REAL (KIND=nag_wp), INTENT (IN) :: s, y ! .. Executable Statements .. cg1 = y RETURN END FUNCTION cg1 FUNCTION ck2(t) ! .. Use Statements .. USE nag_library, ONLY : x01aaf ! .. Implicit None Statement .. IMPLICIT NONE ! .. Function Return Value .. REAL (KIND=nag_wp) :: ck2 ! .. Scalar Arguments .. REAL (KIND=nag_wp), INTENT (IN) :: t ! .. Local Scalars .. REAL (KIND=nag_wp) :: pi ! .. Intrinsic Functions .. INTRINSIC sqrt ! .. Executable Statements .. ck2 = sqrt(x01aaf(pi)) RETURN END FUNCTION ck2 FUNCTION cf2(t) ! .. Implicit None Statement .. IMPLICIT NONE ! .. Function Return Value .. REAL (KIND=nag_wp) :: cf2 ! .. Scalar Arguments .. REAL (KIND=nag_wp), INTENT (IN) :: t ! .. Local Scalars .. REAL (KIND=nag_wp) :: st1 ! .. Intrinsic Functions .. INTRINSIC log, sqrt ! .. Executable Statements .. st1 = sqrt(1.0_nag_wp+t) cf2 = -2.0_nag_wp*log(st1+sqrt(t))/st1 RETURN END FUNCTION cf2 FUNCTION cg2(s,y) ! .. Implicit None Statement .. IMPLICIT NONE ! .. Function Return Value .. REAL (KIND=nag_wp) :: cg2 ! .. Scalar Arguments .. REAL (KIND=nag_wp), INTENT (IN) :: s, y ! .. Executable Statements .. cg2 = y RETURN END FUNCTION cg2 END MODULE d05befe_mod PROGRAM d05befe ! D05BEF Example Main Program ! .. Use Statements .. USE nag_library, ONLY : d05bef, x02ajf USE d05befe_mod, ONLY : cf1, cf2, cg1, cg2, ck1, ck2, iorder, lct, lwk, & nag_wp, nmesh, nout, sol1, sol2 ! .. Implicit None Statement .. IMPLICIT NONE ! .. Local Scalars .. REAL (KIND=nag_wp) :: err, errmax, h, hi1, soln, t, & tlim, tolnl INTEGER :: i, ifail ! .. Local Arrays .. REAL (KIND=nag_wp) :: work(lwk), yn(nmesh) INTEGER :: nct(lct) ! .. Intrinsic Functions .. INTRINSIC abs, mod, real, sqrt ! .. Executable Statements .. WRITE (nout,*) 'D05BEF Example Program Results' tlim = 7.0_nag_wp tolnl = sqrt(x02ajf()) h = tlim/real(nmesh-1,kind=nag_wp) yn(1) = 0.0_nag_wp ifail = 0 CALL d05bef(ck1,cf1,cg1,'Initial',iorder,tlim,tolnl,nmesh,yn,work,lwk, & nct,ifail) WRITE (nout,*) WRITE (nout,*) 'Example 1' WRITE (nout,*) WRITE (nout,99997) h WRITE (nout,*) WRITE (nout,*) ' T Approximate' WRITE (nout,*) ' Solution ' WRITE (nout,*) errmax = 0.0_nag_wp DO i = 2, nmesh hi1 = real(i-1,kind=nag_wp)*h err = abs(yn(i)-sol1(hi1)) IF (err>errmax) THEN errmax = err t = hi1 soln = yn(i) END IF IF (i>5 .AND. mod(i,5)==1) THEN WRITE (nout,99998) hi1, yn(i) END IF END DO WRITE (nout,*) WRITE (nout,99999) errmax, t, soln WRITE (nout,*) tlim = 5.0_nag_wp h = tlim/real(nmesh-1,kind=nag_wp) yn(1) = 1.0_nag_wp ifail = 0 CALL d05bef(ck2,cf2,cg2,'Subsequent',iorder,tlim,tolnl,nmesh,yn,work, & lwk,nct,ifail) WRITE (nout,*) WRITE (nout,*) 'Example 2' WRITE (nout,*) WRITE (nout,99997) h WRITE (nout,*) WRITE (nout,*) ' T Approximate' WRITE (nout,*) ' Solution ' WRITE (nout,*) errmax = 0.0_nag_wp DO i = 1, nmesh hi1 = real(i-1,kind=nag_wp)*h err = abs(yn(i)-sol2(hi1)) IF (err>errmax) THEN errmax = err t = hi1 soln = yn(i) END IF IF (i>7 .AND. mod(i,7)==1) THEN WRITE (nout,99998) hi1, yn(i) END IF END DO WRITE (nout,*) WRITE (nout,99999) errmax, t, soln 99999 FORMAT (' The maximum absolute error, ',E10.2,', occurred at T =', & F8.4/' with solution ',F8.4) 99998 FORMAT (1X,F8.4,F15.4) 99997 FORMAT (' The stepsize h = ',F8.4) END PROGRAM d05befe