! D02MZF Example Program Text ! Mark 24 Release. NAG Copyright 2012. Module d02mzfe_mod ! D02MZF 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 :: tstep = 0.02_nag_wp Integer, Parameter :: iset = 1, neq = 3, nin = 5, & nout = 6 Integer, Parameter :: nrw = 50 + 4*neq Integer, Parameter :: nwkjac = neq*(neq+1) Integer, Parameter :: sdysav = 6 Integer, Parameter :: ldysav = neq Contains Subroutine resid(neq,t,y,ydot,r,ires) ! .. Scalar Arguments .. Real (Kind=nag_wp), Intent (In) :: t Integer, Intent (Inout) :: ires Integer, Intent (In) :: neq ! .. Array Arguments .. Real (Kind=nag_wp), Intent (Out) :: r(neq) Real (Kind=nag_wp), Intent (In) :: y(neq), ydot(neq) ! .. Local Scalars .. Real (Kind=nag_wp) :: alpha = 0.04_nag_wp Real (Kind=nag_wp) :: beta = 1.0E4_nag_wp Real (Kind=nag_wp) :: gamma = 3.0E7_nag_wp ! .. Executable Statements .. r(1:3) = -ydot(1:3) If (ires==1) Then r(1) = -alpha*y(1) + beta*y(2)*y(3) + r(1) r(2) = alpha*y(1) - beta*y(2)*y(3) - gamma*y(2)*y(2) + r(2) r(3) = gamma*y(2)*y(2) + r(3) End If Return End Subroutine resid End Module d02mzfe_mod Program d02mzfe ! D02MZF Example Main Program ! .. Use Statements .. Use nag_library, Only: d02mzf, d02nby, d02ngf, d02ngz, d02nsf, d02nvf, & d02nyf, nag_wp, x04abf Use d02mzfe_mod, Only: iset, ldysav, neq, nin, nout, nrw, nwkjac, resid, & sdysav, tstep ! .. Implicit None Statement .. Implicit None ! .. Local Scalars .. Real (Kind=nag_wp) :: h, h0, hmax, hmin, hu, t, tcrit, & tcur, tolsf, tout, xout Integer :: i, ifail, imxer, iout, itask, & itol, itrace, maxord, maxstp, & mxhnil, niter, nje, nq, nqu, & nre, nst, outchn, pstat Logical :: petzld ! .. Local Arrays .. Real (Kind=nag_wp), Allocatable :: atol(:), rtol(:), rwork(:), & sol(:), wkjac(:), y(:), ydot(:), & ysav(:,:) Real (Kind=nag_wp) :: con(6) Integer :: inform(23) Logical, Allocatable :: algequ(:) Logical :: lderiv(2) ! .. Executable Statements .. Write (nout,*) 'D02MZF Example Program Results' ! Skip heading in data file Read (nin,*) ! Allocations based on number of differential equations (neq) Allocate (atol(neq),rtol(neq),rwork(nrw),sol(neq),wkjac(nwkjac),y(neq), & ydot(neq),ysav(ldysav,sdysav),algequ(neq)) ! Read algorithmic parameters Read (nin,*) maxord, maxstp, mxhnil, pstat Read (nin,*) petzld Read (nin,*) hmin, hmax, h0, tcrit Read (nin,*) t, tout Read (nin,*) itol Read (nin,*) y(1:neq) Read (nin,*) lderiv(1:2) Read (nin,*) rtol(1:neq) Read (nin,*) atol(1:neq) outchn = nout Call x04abf(iset,outchn) con(1:6) = 0.0_nag_wp itask = 2 itrace = 0 ! ifail: behaviour on error exit ! =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft ifail = 0 Call d02nvf(neq,sdysav,maxord,'Functional-iteration',petzld,con,tcrit, & hmin,hmax,h0,maxstp,mxhnil,'Average-L2',rwork,ifail) ! Linear algebra setup. ifail = 0 Call d02nsf(neq,neq,'Numerical',nwkjac,rwork,ifail) Write (nout,99994)(i,i=1,neq) Write (nout,99999) t, y(1:neq) ! First value of t to interpolate and print solution at. iout = 1 xout = tstep ! Integrate one step at a time and overshoot tout (itask=2). steps: Do ifail = 0 Call d02ngf(neq,ldysav,t,tout,y,ydot,rwork,rtol,atol,itol,inform, & resid,ysav,sdysav,d02ngz,wkjac,nwkjac,d02nby,lderiv,itask,itrace, & ifail) ! Get Current value of t (tcur) and last step used (hu) Call d02nyf(neq,neq,hu,h,tcur,tolsf,rwork,nst,nre,nje,nqu,nq,niter, & imxer,algequ,inform,ifail) interp: Do If (tcur-hu=6) Then ! Final solution point printed. Print stats if required. If (pstat==1) Then Write (nout,*) Write (nout,99998) hu, h, tcur Write (nout,99997) nst, nre, nje Write (nout,99996) nqu, nq, niter Write (nout,99995) ' Max err comp = ', imxer End If Exit steps End If ! Set next xout. This might also be in last step, so keep ! looping for interpolation. xout = xout + tstep Cycle interp Else ! Take another step. Cycle steps End If End Do interp End Do steps 99999 Format (1X,F8.3,3(F13.5,2X)) 99998 Format (1X,' HUSED = ',E12.5,' HNEXT = ',E12.5,' TCUR = ',E12.5) 99997 Format (1X,' NST = ',I6,' NRE = ',I6,' NJE = ',I6) 99996 Format (1X,' NQU = ',I6,' NQ = ',I6,' NITER = ',I6) 99995 Format (1X,A,I4) 99994 Format (/1X,' X ',3(' Y(',I1,') ')) End Program d02mzfe