! D05BEF Example Program Text
! Mark 25 Release. NAG Copyright 2014.
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
! .. Accessibility Statements ..
Private
Public :: cf1, cf2, cg1, cg2, ck1, ck2, &
sol1, sol2
! .. Parameters ..
Integer, Parameter, Public :: iorder = 4
Integer, Parameter, Public :: nmesh = 2**6 + 2*iorder - 1
Integer, Parameter, Public :: nout = 6
Integer, Parameter, Public :: lct = nmesh/32 + 1
Integer, Parameter, Public :: &
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