! D03PEF Example Program Text ! Mark 24 Release. NAG Copyright 2012. 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 ! .. 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 Procedures .. 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) ! .. 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) ! .. 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 Procedures .. 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) ! .. 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 Procedures .. 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 Procedures .. 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