! D02MZF Example Program Text
! Mark 25 Release. NAG Copyright 2014.
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
! .. Accessibility Statements ..
Private
Public :: resid
! .. Parameters ..
Real (Kind=nag_wp), Parameter, Public :: tstep = 0.02_nag_wp
Integer, Parameter, Public :: iset = 1, neq = 3, nin = 5, &
nout = 6
Integer, Parameter, Public :: nrw = 50 + 4*neq
Integer, Parameter, Public :: nwkjac = neq*(neq+1)
Integer, Parameter, Public :: sdysav = 6
Integer, Parameter, Public :: ldysav = neq
Contains
Subroutine resid(neq,t,y,ydot,r,ires)
! .. Parameters ..
Real (Kind=nag_wp), Parameter :: alpha = 0.04_nag_wp
Real (Kind=nag_wp), Parameter :: beta = 1.0E4_nag_wp
Real (Kind=nag_wp), Parameter :: gamma = 3.0E7_nag_wp
! .. 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)
! .. 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<xout .And. xout<=tcur) Then
! xout is in interval of last step, so interpolate.
ifail = 0
Call d02mzf(xout,sol,neq,ldysav,neq,ysav,sdysav,rwork,ifail)
Write (nout,99999) xout, sol(1:neq)
iout = iout + 1
If (iout>=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