! D05BEF Example Program Text ! Mark 24 Release. NAG Copyright 2012. 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 ! .. 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 Procedures .. 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) ! .. 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) ! .. Function Return Value .. Real (Kind=nag_wp) :: ck1 ! .. Scalar Arguments .. Real (Kind=nag_wp), Intent (In) :: t ! .. Intrinsic Procedures .. 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 ! .. 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 Procedures .. 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) ! .. 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 ! .. 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 Procedures .. Intrinsic :: sqrt ! .. Executable Statements .. ck2 = sqrt(x01aaf(pi)) Return End Function ck2 Function cf2(t) ! .. 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 Procedures .. 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) ! .. 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, nag_wp, x02ajf Use d05befe_mod, Only: cf1, cf2, cg1, cg2, ck1, ck2, iorder, lct, lwk, & 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 Procedures .. 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