! D02ZAF Example Program Text
! Mark 26.1 Release. NAG Copyright 2016.
Module d02zafe_mod
! D02ZAF 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 :: fcn, jac, monitr
! .. Parameters ..
Integer, Parameter, Public :: iset = 1, itrace = 0, neq = 3, &
nin = 5, nout = 6
Integer, Parameter, Public :: nrw = 50 + 4*neq
Integer, Parameter, Public :: nwkjac = neq*(neq+1)
Integer, Parameter, Public :: ldysav = neq
Contains
Subroutine fcn(neq,t,y,f,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) :: f(neq)
Real (Kind=nag_wp), Intent (In) :: y(neq)
! .. Executable Statements ..
f(1) = -0.04E0_nag_wp*y(1) + 1.0E4_nag_wp*y(2)*y(3)
f(2) = 0.04E0_nag_wp*y(1) - 1.0E4_nag_wp*y(2)*y(3) - &
3.0E7_nag_wp*y(2)*y(2)
f(3) = 3.0E7_nag_wp*y(2)*y(2)
Return
End Subroutine fcn
Subroutine jac(neq,t,y,h,d,p)
! .. Scalar Arguments ..
Real (Kind=nag_wp), Intent (In) :: d, h, t
Integer, Intent (In) :: neq
! .. Array Arguments ..
Real (Kind=nag_wp), Intent (Inout) :: p(neq,neq)
Real (Kind=nag_wp), Intent (In) :: y(neq)
! .. Local Scalars ..
Real (Kind=nag_wp) :: hxd
! .. Executable Statements ..
hxd = h*d
p(1,1) = 1.0E0_nag_wp - hxd*(-0.04E0_nag_wp)
p(1,2) = -hxd*(1.0E4_nag_wp*y(3))
p(1,3) = -hxd*(1.0E4_nag_wp*y(2))
p(2,1) = -hxd*(0.04E0_nag_wp)
p(2,2) = 1.0E0_nag_wp - hxd*(-1.0E4_nag_wp*y(3)-6.0E7_nag_wp*y(2))
p(2,3) = -hxd*(-1.0E4_nag_wp*y(2))
! Do not need to set P(3,1) since Jacobian preset to zero
! P(3,1) = - HXD*(0.0E0)
p(3,2) = -hxd*(6.0E7_nag_wp*y(2))
p(3,3) = 1.0E0_nag_wp - hxd*(0.0E0_nag_wp)
Return
End Subroutine jac
Subroutine monitr(neq,ldysav,t,hlast,hnext,y,ydot,ysav,r,acor,imon,inln, &
hmin,hmax,nqu)
! .. Use Statements ..
Use nag_library, Only: d02zaf
! .. Scalar Arguments ..
Real (Kind=nag_wp), Intent (In) :: hlast, t
Real (Kind=nag_wp), Intent (Inout) :: hmax, hmin, hnext
Integer, Intent (Inout) :: imon
Integer, Intent (Out) :: inln
Integer, Intent (In) :: ldysav, neq, nqu
! .. Array Arguments ..
Real (Kind=nag_wp), Intent (In) :: acor(neq,2), r(neq), ydot(neq), &
ysav(ldysav,*)
Real (Kind=nag_wp), Intent (Inout) :: y(neq)
! .. Local Scalars ..
Real (Kind=nag_wp) :: errloc
Integer :: i, ifail
! .. Executable Statements ..
inln = 3
If (imon==1) Then
ifail = -1
errloc = d02zaf(neq,acor(1,2),acor(1,1),ifail)
If (ifail/=0) Then
imon = -2
Else If (errloc>5.0E0_nag_wp) Then
Write (nout,99999) t, (y(i),i=1,neq), errloc
Else
Write (nout,99998) t, (y(i),i=1,neq)
End If
End If
Return
99999 Format (1X,F10.6,3(F13.7,2X),/,1X,' ** WARNING scaled local error = ', &
F13.5)
99998 Format (1X,F10.6,3(F13.7,2X))
End Subroutine monitr
End Module d02zafe_mod
Program d02zafe
! D02ZAF Example Main Program
! .. Use Statements ..
Use d02zafe_mod, Only: fcn, iset, itrace, jac, ldysav, monitr, neq, nin, &
nout, nrw, nwkjac
Use nag_library, Only: d02nbf, d02nsf, d02nvf, d02nyf, nag_wp, x04abf
! .. Implicit None Statement ..
Implicit None
! .. Local Scalars ..
Real (Kind=nag_wp) :: h, h0, hmax, hmin, hu, t, tcrit, &
tcur, tolsf, tout
Integer :: i, ifail, imxer, itask, itol, &
maxord, maxstp, mxhnil, niter, nje, &
nq, nqu, nre, nst, outchn, sdysav
Logical :: petzld
! .. Local Arrays ..
Real (Kind=nag_wp), Allocatable :: atol(:), rtol(:), rwork(:), &
wkjac(:), y(:), ydot(:), ysav(:,:)
Real (Kind=nag_wp) :: con(6)
Integer :: inform(23)
Logical, Allocatable :: algequ(:)
! .. Executable Statements ..
Write (nout,*) 'D02ZAF Example Program Results'
! Skip heading in data file
Read (nin,*)
! neq: number of differential equations
Read (nin,*) maxord, maxstp, mxhnil
sdysav = maxord + 1
Allocate (atol(neq),rtol(neq),rwork(nrw),wkjac(nwkjac),y(neq),ydot(neq), &
ysav(ldysav,sdysav),algequ(neq))
outchn = nout
Call x04abf(iset,outchn)
! Set algorithmic and problem parameters
Read (nin,*) hmin, hmax, h0, t, tout
Read (nin,*) petzld
! Initialization
! Integrate to tout without passing tout.
tcrit = tout
itask = 4
! Use default values for the array con.
con(1:6) = 0.0_nag_wp
! Use BDF formulae with modified Newton method.
! Use averaged L2 norm for local error control.
ifail = 0
Call d02nvf(neq,sdysav,maxord,'Newton',petzld,con,tcrit,hmin,hmax,h0, &
maxstp,mxhnil,'Average-L2',rwork,ifail)
! Setup for using analytic Jacobian
ifail = 0
Call d02nsf(neq,neq,'Analytical',nwkjac,rwork,ifail)
Write (nout,*)
Write (nout,*) ' Analytic Jacobian'
Write (nout,*)
! Set tolerances.
Read (nin,*) itol
Read (nin,*) rtol(1), atol(1)
! Initial values for Y.
Read (nin,*) y(1:neq)
Write (nout,*) ' X Y(1) Y(2) Y(3)'
Write (nout,99999) t, (y(i),i=1,neq)
! Solve the problem using MONITR to output intermediate results.
ifail = -1
Call d02nbf(neq,ldysav,t,tout,y,ydot,rwork,rtol,atol,itol,inform,fcn, &
ysav,sdysav,jac,wkjac,nwkjac,monitr,itask,itrace,ifail)
If (ifail==0) Then
! Get integration statistics.
Call d02nyf(neq,neq,hu,h,tcur,tolsf,rwork,nst,nre,nje,nqu,nq,niter, &
imxer,algequ,inform,ifail)
Write (nout,*)
Write (nout,99997) hu, h, tcur
Write (nout,99996) nst, nre, nje
Write (nout,99995) nqu, nq, niter
Write (nout,99994) ' Max Err Comp = ', imxer
Write (nout,*)
Else
Write (nout,*)
Write (nout,99998) 'Exit D02NBF with IFAIL = ', ifail, ' and T = ', t
End If
99999 Format (1X,F10.6,3(F13.7,2X))
99998 Format (1X,A,I2,A,E12.5)
99997 Format (1X,' HUSED = ',E12.5,' HNEXT = ',E12.5,' TCUR = ',E12.5)
99996 Format (1X,' NST = ',I6,' NRE = ',I6,' NJE = ',I6)
99995 Format (1X,' NQU = ',I6,' NQ = ',I6,' NITER = ',I6)
99994 Format (1X,A,I4)
End Program d02zafe