! D03PEF Example Program Text
! Mark 23 Release. NAG Copyright 2011.
MODULE d03pefe_mod
! D03PEF Example Program Module:
! Parameters and User-defined Routines
! .. Use Statements ..
USE nag_library, ONLY : nag_wp
! .. Implicit None Statement ..
IMPLICIT NONE
! .. Parameters ..
REAL (KIND=nag_wp), PARAMETER :: half = 0.5_nag_wp
REAL (KIND=nag_wp), PARAMETER :: one = 1.0_nag_wp
INTEGER, PARAMETER :: nin = 5, nleft = 1, nout = 6, &
npde = 2
CONTAINS
SUBROUTINE uinit(npde,npts,x,u)
! Routine for PDE initial values
! .. Implicit None Statement ..
IMPLICIT NONE
! .. Scalar Arguments ..
INTEGER, INTENT (IN) :: npde, npts
! .. Array Arguments ..
REAL (KIND=nag_wp), INTENT (OUT) :: u(npde,npts)
REAL (KIND=nag_wp), INTENT (IN) :: x(npts)
! .. Local Scalars ..
INTEGER :: i
! .. Intrinsic Functions ..
INTRINSIC exp, sin
! .. Executable Statements ..
DO i = 1, npts
u(1,i) = exp(x(i))
u(2,i) = sin(x(i))
END DO
RETURN
END SUBROUTINE uinit
SUBROUTINE pdedef(npde,t,x,u,ut,ux,res,ires)
! .. Implicit None Statement ..
IMPLICIT NONE
! .. Scalar Arguments ..
REAL (KIND=nag_wp), INTENT (IN) :: t, x
INTEGER, INTENT (INOUT) :: ires
INTEGER, INTENT (IN) :: npde
! .. Array Arguments ..
REAL (KIND=nag_wp), INTENT (OUT) :: res(npde)
REAL (KIND=nag_wp), INTENT (IN) :: u(npde), ut(npde), ux(npde)
! .. Executable Statements ..
IF (ires==-1) THEN
res(1) = ut(1)
res(2) = ut(2)
ELSE
res(1) = ut(1) + ux(1) + ux(2)
res(2) = ut(2) + 4.0_nag_wp*ux(1) + ux(2)
END IF
RETURN
END SUBROUTINE pdedef
SUBROUTINE bndary(npde,t,ibnd,nobc,u,ut,res,ires)
! .. Implicit None Statement ..
IMPLICIT NONE
! .. Scalar Arguments ..
REAL (KIND=nag_wp), INTENT (IN) :: t
INTEGER, INTENT (IN) :: ibnd, nobc, npde
INTEGER, INTENT (INOUT) :: ires
! .. Array Arguments ..
REAL (KIND=nag_wp), INTENT (OUT) :: res(nobc)
REAL (KIND=nag_wp), INTENT (IN) :: u(npde), ut(npde)
! .. Local Scalars ..
REAL (KIND=nag_wp) :: t1, t3
! .. Intrinsic Functions ..
INTRINSIC exp, sin
! .. Executable Statements ..
IF (ires==-1) THEN
res(1) = 0.0_nag_wp
ELSE IF (ibnd==0) THEN
t3 = -3.0_nag_wp*t
t1 = t
res(1) = u(1) - half*((exp(t3)+exp(t1))+half*(sin(t3)-sin(t1)))
ELSE
t3 = one - 3.0_nag_wp*t
t1 = one + t
res(1) = u(2) - ((exp(t3)-exp(t1))+half*(sin(t3)+sin(t1)))
END IF
RETURN
END SUBROUTINE bndary
SUBROUTINE exact(t,npde,npts,x,u)
! Exact solution (for comparison purposes)
! .. Implicit None Statement ..
IMPLICIT NONE
! .. Scalar Arguments ..
REAL (KIND=nag_wp), INTENT (IN) :: t
INTEGER, INTENT (IN) :: npde, npts
! .. Array Arguments ..
REAL (KIND=nag_wp), INTENT (OUT) :: u(npde,npts)
REAL (KIND=nag_wp), INTENT (IN) :: x(npts)
! .. Local Scalars ..
REAL (KIND=nag_wp) :: xt, xt3
INTEGER :: i
! .. Intrinsic Functions ..
INTRINSIC exp, sin
! .. Executable Statements ..
DO i = 1, npts
xt3 = x(i) - 3.0_nag_wp*t
xt = x(i) + t
u(1,i) = half*((exp(xt3)+exp(xt))+half*(sin(xt3)-sin(xt)))
u(2,i) = (exp(xt3)-exp(xt)) + half*(sin(xt3)+sin(xt))
END DO
RETURN
END SUBROUTINE exact
END MODULE d03pefe_mod
PROGRAM d03pefe
! D03PEF Example Main Program
! .. Use Statements ..
USE nag_library, ONLY : d03pef, nag_wp
USE d03pefe_mod, ONLY : bndary, exact, nin, nleft, nout, npde, pdedef, &
uinit
! .. Implicit None Statement ..
IMPLICIT NONE
! .. Local Scalars ..
REAL (KIND=nag_wp) :: acc, tout, ts
INTEGER :: i, ifail, ind, it, itask, &
itrace, lisave, lrsave, neqn, &
npts, nwkres
! .. Local Arrays ..
REAL (KIND=nag_wp), ALLOCATABLE :: eu(:,:), rsave(:), u(:,:), x(:)
INTEGER, ALLOCATABLE :: isave(:)
! .. Intrinsic Functions ..
INTRINSIC real
! .. Executable Statements ..
WRITE (nout,*) 'D03PEF Example Program Results'
! Skip heading in data file
READ (nin,*)
READ (nin,*) npts
nwkres = npde*(npts+21+3*npde) + 7*npts + 4
neqn = npde*npts
lisave = neqn + 24
lrsave = 11*neqn + (4*npde+nleft+2)*neqn + 50 + nwkres
ALLOCATE (eu(npde,npts),rsave(lrsave),u(npde,npts),x(npts), &
isave(lisave))
READ (nin,*) acc
READ (nin,*) itrace
! Set spatial-mesh points
DO i = 1, npts
x(i) = real(i-1,kind=nag_wp)/real(npts-1,kind=nag_wp)
END DO
ind = 0
itask = 1
CALL uinit(npde,npts,x,u)
! Loop over output value of t
READ (nin,*) ts, tout
DO it = 1, 5
tout = 0.2_nag_wp*real(it,kind=nag_wp)
! ifail: behaviour on error exit
! =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft
ifail = 0
CALL d03pef(npde,ts,tout,pdedef,bndary,u,npts,x,nleft,acc,rsave, &
lrsave,isave,lisave,itask,itrace,ind,ifail)
IF (it==1) THEN
WRITE (nout,99997) acc, npts
WRITE (nout,99999) x(5), x(13), x(21), x(29), x(37)
END IF
! Check against the exact solution
CALL exact(tout,npde,npts,x,eu)
WRITE (nout,99998) ts
WRITE (nout,99995) u(1,5:37:8)
WRITE (nout,99994) eu(1,5:37:8)
WRITE (nout,99993) u(2,5:37:8)
WRITE (nout,99992) eu(2,5:37:8)
END DO
WRITE (nout,99996) isave(1), isave(2), isave(3), isave(5)
99999 FORMAT (' X ',5F10.4/)
99998 FORMAT (' T = ',F5.2)
99997 FORMAT (//' Accuracy requirement =',E10.3,' Number of points = ',I3/)
99996 FORMAT (' Number of integration steps in time = ',I6/' Number o', &
'f function evaluations = ',I6/' Number of Jacobian eval', &
'uations =',I6/' Number of iterations = ',I6)
99995 FORMAT (' Approx U1',5F10.4)
99994 FORMAT (' Exact U1',5F10.4)
99993 FORMAT (' Approx U2',5F10.4)
99992 FORMAT (' Exact U2',5F10.4/)
END PROGRAM d03pefe