! 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